NBS-GCR-83-451 


Computer Fire Code VI - Volume 1 


_ December 1983 


Sponsored by 

U.S. DEPARTMENT OF COMMERCE 
National Bureau of Standards 

Center for Fire Research 

Washington, DC 20234 


Aw abr e 
a ee * : ik 
a ae 7 ; Ti) 
t io ve : a Al 
det en 7 
* . ; 
b = ae | 
‘ 
z 
oP 
7 
" 
af = ; 
7 al 
i 7 a 
: ig ~ 
A! | ; “a 
; To _ oe al 


y 
ss oe 
; i 
Ss Med “15 
a ry 


ry! 


a 
me 
‘ 


NBS-GCR-83-451 


COMPUTER FIRE CODE VI - VOLUME 1 


J.B. Gahm 


Harvard University 
Division of Applied Sciences 
Cambridge, MA 


December 1983 


Sponsored by 

U.S. DEPARTMENT OF COMMERCE 
National Bureau of Standards 
Center for Fire Research 
Washington, DC 20234 


NOTICE 


This report was prepared for the Center for Fire 
Research of the National Engineering Laboratory, 
National Bureau of Standards under Grant Number 
NB82NADA3030. The statements and conclusions 
contained in this report are those of the authors 
and do not necessarily reflect the views of the 
National Bureau of Standards or the Center for 
Fire Research. , 


Li 


Computer Fire Code VI 


J. B. Gahm 


Division of Applied Sciences 
July 1983 


Home Fire Project Technical Report No. 58 


iii 


Table of Contents 


Abstract’ . .¢ «Nab ele ie es 8S ee eae eke been, ce eee 
FODEWOLOs. tbe oo Pe ear Oe Re 
Chapter 1. Introduction’. 2... 00505 6 0st et rte eo eee 1 
1.1 Intended Audience” °°.’.”, ¢#@f vice oh. SORE o",.. 1 
1.2 Overview of the Structure’ 23°) 7802, Ge SQhc hasty. ee 2 
1.2.1 Physical Equations and Physical Subroutines. ......... 2 
1.2.2 The System of Equations. ww %) cnc: «6 sue wi nun oni lei ee 3 
1.2.8 The,Numerics Module <p. cc ioc oss one, meeps isnen nn 5 
1.2.4 The Interface Module .... 2.540 5 0+ & eisseis einen 5 
1.2.4.1 The Function of the Interface. ..........2000- 5 
1.2.4.2 The Structure of the Interface ...........2000- T 
Chapter 2. The Physics Module....... ie vee a 8 6 a ae 8 
2.1 The Physical Equations ....... si... 5 « «0s. see 8 
2.1.1 A Single Equation in the System of Equations......... 8 
2.1.2 The Forms of a Physical Equation ..........6.24-. 11 
2.1.3 Restrictions on the Physical Equations ............ 13 
2.2 The Physical Subroutines 2 9.02 2)..0¢. 6) co.cc eee 15 
2.2.1 Physical Subroutines and Physical Equations ........ 15 
2.2.2 The Structure Of a Physical Subroutine ........... 17 
2:2.2.1. “inputs and Outputs =... ee 17 
2.2.2.2 Data Structures Within a Physical Subroutine. .... 18 
2.2.2.3 Coding Standards for a Physical Subroutine ...... 19 
2.2.2.4 Non-Functional Subroutines ............206-. 21 
2.3 Thoughts for Future Development .............-008- 24 
Chapter 3: The: Numerics Module .. 20.5. « . 0.3) 5 eee ee 25 
3.1 ‘Subroutine ALDIFF os 0. 2.5 oe oe 25 
3.1.1 Starting the Integration 4.2.5...) ae ee 26 
3.1.2 A Normal Integration Cycle ....... rn te 26 
3.1.3 Notes on the Internal Workings of ALDIFF........... 28 
3.2 Restrictions on the System of Equations .............. 29 


iv 


Table of Contents 


henier saebnenitertace:MOdUIC. «fi, oo sues. goss vate ees. vie eowles 30 
4.1 The Control Section of the Interface. ............... 32 

4 leverProgram MARKG: wotewev. avi nieltes won . GOORL nade) us 32 

4.12. SUDIOULING UNPINT o,f FOL O26P yee YO? Bitay of Lifes, Jay 36 

A ASMPOUDLOULINCMMMS TAL cb: ge cies chs con oid. nos. t auto sath ie 39 

41:42 Subroutine SETSYS harsses vlearet. « ob. a0 st bawheters aa ad 44 

eres MOUDrOULINe NENSTP eo%*, At wetects .eds, Al. JonvedaAl adlaye 48 

4.2 The Interface Communication Section............... 49 
AsD bei PDEGUiINC URBACK ti. as ectt 's aienduey canes wash? east CTS 49 

2.28 Subroutine FURCTS Ia! 02 wow, 22. Jade. earl. bagnedoay.et 29: 49 

4.2.3 Subroutines EVALFP and ERFDE............-e--ce0- 51 

Ae2 Say ariaviesGMinaulOl 2c erly pes eld sos 2 nese 008 51 

4.2.3.2 Organization of the Subroutines ............. 56 
Appendix A: Introducing a New Physical Subroutine ........... 58 
Modifying the Existing Equations. ........0.-¢,7.00+s808085 58 
Adding a New PhysicalSubroutine shies. au 0 siemens cer erecey ole pose 59 
Oia Newer ivcical Variable . oc. . csc sis 2s eos 9s peierre oper 59 
Appendix B: Numerical Considerations for Fire Models ......... 62 
LDL TOUMGM OTe MEN PSta%G. 64 . 1.90) sededad-mowt> nosed 9%) na9 62 
CEO OTe as te AS te ARE POT AR LTR or Tega etek 64 
NE Ne sca tca a ous tt ene tees 66 

Pe ILCENALIMGOOMILION MCLDOG! 2. oss teoe,> ciceta'sdecs tend stiemen eh ois 67 
PeRPENOIS 7. LICHIOUET Ys Oe 2%, 020) oo tPRY, Ot ott wale) te) tietow wor ook 68 
Sedat EY STE gta Gobel Manali tei mete baer tn tatarll tut erie belle OO le ee he 104 
Pr CrcreliCeSanhes hed Seed (alse rds aikow shisresals +. N+ VES taste) tantatutart 105 


Appendix D appears in a separate volume: 
Part 1 - BINPS 
Part 2 - DBLES 
Part 3 - Common Blocks 


ABSTRACT 


This is the documentation for computer Fire Code VI (CFC VI), which 
was originally meant to be the successor to CFC V. There is, as of this 
writing (November 1983), no definitive version of CFC VI. However, this 
document will be valid for any version, since it mostly describes the 
structure of the program, which will not change. The numerical "package" 
used to solve the equations is described in general terms; so far as 
details are concerned, it is largely treated as if it were a "black box." 
The physics inherent in the program is almost entirely omitted; however, 
most of the subroutines giving the physics have been taken directly from 
Mark 5, although they have been broken up so that any subroutine (or 
FUNCTION) has just one output. Hence, much of the documentation in the 
listing is unchanged from what it was in Mark 5. For the reader interest- 
ed in the documentation of the physics, it is given in detail in Harvard 
Fire Project Technical Reports 34 and 45. 


There are three important appendices, wherein it is explained, for 
anyone who may wish to develop the program further, how to insert a new 
subroutine or a new physical variable; also, a dictionary of symbols, 
terms, and variables used in the program and in the text. A few 
additional germane comments are made in the Foreword. 


Following this report, a listing of a version of CFC VI appears; 
this is a version which compiles on both a VAX and a Perkin-Elmer 
computer (it dates from October 1983). The program was written in 
(approximately) ANSI 77 FORTRAN, and is fairly machine-independent. 
The program consists of two independent parts: BINP is a program which 
produces an input file to be used by the second (main) program, DBLE. 


From the point of view of the user, CFC VI differs from Mark 5 
in a number of respects, all described in the Foreword except for 
two points: first, the gas burner algorithm in Mark 5 has not been 
incorporated into CFC VI. Second, some error(s) have been made in 
transcribing the species concentratin algorithm from 5 to 6, so that 
the gas concentrations (and perhaps that of smoke as well) are in- 
correct. However, the error(s) should be easy to find and correct. 
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FOREWORD 


by H. E. Mitler 


This documentation of CFC VI is independent of that for CFC V, except 


that the physical equations of the model are described in the earlier document, 
whereas they are not described here. Reading the earlier documentation might 
also make reading this a bit easier. 


CFC VI incorporates features which make it superior to CFC V in certain 


ways: 


Le 


The structure is more modular, with the physics and numerics more nearly 
separated. 


2. A scheme for “elimination” of dependent variables has been (implicitly) 


incorporated, so that the principal numerical manipulation occurs for a much 
reduced set of equations. 


It has been generalized into a multi-room (on one floor) program. 


The numerical package guarantees convergence within a predetermined ac- 
curacy (so long as the equations satisfy certain criteria). With the existing 
equations, it is (in some ways) in fact more robust, so that it often converges 
to a solution at points where CFC V fails to converge. 


The input section has been split off from the main program, and made 
more “user-friendly:” it gives “random access” to input data that needs 
to be changed, rather than having to go through all of the default values 
serially. It has a built-in handbook of needed thermophysical data. And it 
creates separate input files, any one of which can then be called by the main 
program.! 


1The input package described here has actually been implemented on a version of CFC VI 
which differs slightly from that described in the remainder of this document.—J. Gahm. 


vil 
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There is, however, a certain price to be paid for these advantages: 


First, it takes from three to ten times as long (in CPU time) to make a 
run. There are several reasons for this: the principal one (perhaps) being that 
checking for accuracy is very time-consuming; that sometimes very tiny step 
sizes have to be chosen; and that the program must be run in double precision. 
The contribution made by the different logical flow cannot be estimated. 


Second, sometimes CFC VI fails to converge where CFC V_ has no difficulty. 


Third: although the program was designed top-down, and the logical flow 
should be more transparent, it has in fact turned out in practice to be a more 
difficult program to follow and modify—partly because there has not yet been 
documentation for it, containing (for example) a dictionary of symbols; partly 
because the variable-elimination scheme forces one to call functions and sub- 
routines in a particular order; and partly for other reasons. 


Notwithstanding these drawbacks, the fact that this version of the numerics 
guarantees convergence (when certain restrictions on the physical equations are 
satisfied) and a certain accuracy (per step, of course) is very useful, for it permits 
the program to be used as a reliable standard code. If the listed drawbacks can 
be ameliorated or eliminated, then this could become a production code, as well. 
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Chapter 1 
Introduction 


§1.1 Intended Audience 


This document is intended as programmer rather than user documentation. 
Since CFC VI is intended to be a development code rather than a practical code, 
the users of the program are expected to be somewhat familiar with its basic 
principles of operation. Whereas user documentation would treat the internal 
workings of the program as a “black box,” this document explains what is going 
on inside the box. Conversely, this document omits explanations which would 
seem superfluous to a programmer with the source code in hand. 


Even for the programmer, however, this document cannot stand alone as a 
guide to understanding all the important concepts of CFC VI. First, it makes 
no attempt to explain the physical basis for any of the physical equations which 
comprise the “meat” of the model. In fact, it does not even list the particular 
equations involved in the current model. The physical equations used in CFC VI 
are for the most part the same as those used in CFC V and are given in [4]. 


Second, the discussion of the numerical package (used to solve the fire equa- 
tions) is incomplete. CFC VI uses C. W. Gear’s algebraic and differential 
equation solving algorithm, described in [3]. The discussion of the numerical 
package in this document assumes a general familiarity with Gear’s algorithm 
- and describes only those aspects of the numerical package which are peculiar 
to the CFC VI implementation. Several other discussions in this document also 
assume a familiarity with Gear’s package. 


Third, CFC VI uses a “variable elimination” scheme which is not adequately 
described in this document. “Variable elimination” refers to a procedure which 
substantially reduces the number of equations that the numerical package must 
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solve simultaneously and therefore reduces the amount of computer resources 
needed to run the model. Since much of the structure of CFC VI is designed 
to facilitate this variable elimination scheme, an understanding of that scheme 
should also be considered a prerequisite to understanding CFC VI. Variable 
elimination is fully explained in [1], [6], and [5]. 


Finally, extensive comments included with the source code of CFC VI serve 
as a necessary guide to the details of operation of the program. This document 
outlines only the broader concepts of organization and operation of the program. 


§1.2 Overview of the Structure 


This section attempts to give an overview of the highest level of the structure 
of CFC VI, represented in figure 1. As suggested in the figure, CFC VI is 
conceptually divided into three structural/functional modules, the physics, the 
numerics, and the interface. The physics module embodies a set of equations 
which represent the model’s “knowledge” of the physical processes occurring in 
a fire. The numerics module consists of an integrator which solves systems of 
simultaneous equations from the physics module in order to model a particular 
fire. The interface module is divided into two sub-modules, the control/system 
set-up section and the communications section. The first of these controls 
the integrator and selects from the physics module the appropriate system of 
equations to model a particular fire. The second provides communications 
between the numerics and the physics. The remainder of this section briefly 
describes each of these components of the program and attempts to clarify their 
interrelationships. 


1.2.1 Physical Equations and Physical Subroutines 


The computer model of CFC VI begins with a mathematical model consisting 
of a set of equations which describe the processes being studied. The equations 
of CFC VI consist of algebraic equations and first-order, ordinary differential 
equations in time.! These equations are called physical equations since they 
describe the physical processes of the fire. Each of these physical equations 
is coded in some FORTRAN subroutine called a physical subroutine.? The 


1 at least for the most part. There is also one partial-differential equation, but it requires 
special treatment. See section 2.2.2.4. 


2The same physical subroutine may code more than one physical equation. 


1.2. Overview of the Structure 8 


Saou 


Numerics Module Interface Module Physics Module 


| System 
Set-Up 


Control 


Physical 
Subroutines 


Integrator 


Communications 


Figure 1. Overall Structure of CFC VIL 


collection of all these physical subroutines is called the physics module and is 
represented in the right-hand section of figure 1. Chapter 1 provides a more 
extensive discussion of physical equations, physical subroutines, and how the 
latter implements the former. 


1.2.2 ‘The System of Equations 


Suppose one wants to use CFC VI to model a particular fire. The set of 
physical equations coded in the physics module models a large class of different 
fires. In order to model a particular fire, it is necessary to (1) select those 
equations which are relevant to that fire and (2) specify values for any parameters 
on which those equations might depend. For example, suppose one wishes to 
model a fire involving just a single object. First, one should omit all those 
equations which pertain to other objects. Second, since the physical equations 
are meant to apply to objects with different shapes and physical properties, one 
must specify the particular geometric and physical parameters for this object. 


Selecting the appropriate equations and parameter values for a particular fire 
creates a system of equations corresponding to that fire. A system consists of n 
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simultaneous algebraic and first-order, ordinary differential equations in n vari- 
ables. The program's task is finding the solution to this system of equations— 
finding the value of each of the variables as a function of time. 


The physicist might represent the system of equations simply as a list of 
n equations, using names for the variables and parameters which suggest their 
physical significance. However, representing the equations in a more abstract 
form greatly simplifies the construction of the program. It is worth pausing here 
to look at that form since the functions of both the numerics module and the 
interface module depend on it. 


CFC VI represents the equations of the system in the form 


0 = F(U,V,W,t) (1.1) 

dV /dt = G(U,V,W,t) (1.2) 
W = H(U,V,W,t) (1.3) 

Y = 1(U,V,W,t), (1.4) 


where U, V, W, and Y are vectors of variables; ¢ is the time since the start of 
run; dimn Ff’ = dimnU, dimnG = dimnV, dimnH = dimnW, and dimn/ = 
dimn Y; and dimn F + dimn G+ dimn H +dimn/ = n.° (Here, dimn X means 
the number of components in X.) The individual equations represented by the 
components of (1.2) are the differential equations. The corresponding variables 
(the components of V) are the integrated variables. Those equations represented 
by (1.1), (1.3), and (1.4) are the algebraic equations, and the corresponding 
variables (the components of U, W, and Y) are the algebraic variables. Of 
the algebraic equations, those represented by (1.1) are called root-finding equa- 
tions because any solution of the sytem is a zero or “root” of such an equation. 
The variables corresponding to root-finding equations (the components of U) are 
called root-finding variables. Those equations represented by (1.3) are called 
fixed-point equations because any solution to the system includes a vector, W, 
which appears again (alone) on the left hand side of (1.3). The variables cor- 
responding to fixed-point equations (the components of W) are called fixed-point 
variables. Finally, those equations represented by (1.4) are called eliminated 
variable equations and the corresponding variables (the components of Y) are 
called eliminated variables. Eliminated variables are so called because they have 
been eliminated from the right-hand side of the system of equations. 


3Real physical equations do not depend explicitly on t. Sometimes, the model’s equations 
depend on ? only as a crude approximation to reality. 


Any system of algebraic and first-order, ordinary differential equations can be 
represented in the forms of (1.1)-(1.4). In fact, the first two of these forms would 
suffice, since (1.3) and (1.4) can be readily converted to the form of (1.1). The 
two alternative forms for expressing algebraic equations help to save computation 
time. Any algebraic equation can be easily written as a root-finding equation, but 
CFC VI solves fixed-point equations with less expense and eliminated variable 
equations with almost no expense.* Therefore, one design goal for CFC VI has 
been to find ways to construct the system of equations so as to include as few non- 
eliminated variables as possible. The procedure by which this is accomplished 
is called variable elimination and is discussed in section 4.2.3.1 and [1], [6], and 


[5]. 


1.2.3 The Numerics Module 


The numerics module (represented in the left-most portion of figure 1) solves 
the system of equations, represented in the form of (1.1)-(1.4). It uses an 
algebraic- and differential-equation solver developed by C. W. Gear [3]. Gear’s 
package is a variable order, variable time-step integrator which chooses the order 
of the method and the step size based on estimates of the local truncation error. 
Certain modifications have been made in Gear’s basic package so that it accepts 
the system of equations represented in the form of (1.1)}{1.4). 


For the most part, this document treats Gear’s package as a “black box,” 
relying on Gear’s own documentation to acquaint the reader with its theory and 
structure. Chapter 2 explains the peculiarities of the CFC VI implementation 
of the package. 


1.2.4 The Interface Module 


1.2.4.1 The Function of the Interface 


. The interface must perform three functions: control, system set-up, and 
communications. These are outlined below. 


£7 solve a system which includes many root-finding equations, CFC VI must perform the 
costly operations of evaluating and decomposing the Jacobian matrix of the system. A system 
consisting mostly of fixed-point equations can sometimes be solved by a successive substitution 
method, which avoids these costly operations. Eliminated variables are still better since they 
can be removed entirely from the system of equations to be solved simultaneously: one need 
only solve (1.1), (1.2), and (1.3) simultaneously, substituting the resulting values for U, V, 
and W back in (1.4) to find Y. 
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First, the interface acts as a controller for the integrator. It starts the 
integrator at the beginning of the run, and stops it at the end of the run. 
During the run, it forces the integrator to converge at every time where values 
for the variables are needed. For example, if the user has requested output 
every 10 seconds of fire time, the interface forces the integrator to converge at 
t = 10.0, 20.0, 30.0, ..., etc. 


The second and third functions of the interface require a longer explana- 
tion. Since the system of equations is entirely abstract when represented in 
the forms of (1.1)-(1.4), one must establish some correspondence between the 
real physical equations coded in the physical subroutines and the equations of 
(1.1)-(1.4). For the purposes of the program, this amounts to establishing a cor- 
respondence between the data structures of the numerics and those of the physics. 
The numerics stores values for the variables (the components of (U,V,W,Y)) 
in a one-dimensional array, X, and values for the functions (the components 
of (F,G,H,J)) in a one-dimensional array, Z. It does not have access to any 
parameter values. The physical subroutines, on the other hand, take individual 
variables as input and return individual function values as output. Moreover, 
most require some parameters as input in addition to variables. 


For example, one physical subroutine, CNVL1, requires seven input variables, 
including ZHLZZ (the depth of the hot layer in a room) and ZELZZ (the total 
energy stored in the hot layer of the room), and two input parameters, ZLRZX: 
and ZLRZY (the length and width of the room). It returns one function value 
as output, TELZD1 (the rate of change of layer energy due to convective heating 
of the walls). In order to give physical meaning to the system (1.1)-(1.4), the 
program must somehow interpret between these two different data structures. In 
this case, the program must: (1) decide which positions in X shall contain values 
for ZHLZZ and ZELZZ; (2) decide which (if any) position of Z should receive the 
value output as TELZD1; and (3) supply the parameter values ZLRZX and ZLRZY 
to the physical subroutines. 


The second function of the interface, the system set-up function, is to es- 
tablish this correspondence between the data structures of the numerics and 
those of the physics. As suggested in the preceding example, this involves three 
tasks: (1) decide which position in X corresponds to each of the physical variables; 
(2) decide which (if any) position in Z should receive the output of each physi- 
cal subroutine; and (3) decide which parameter values to pass to each physical 
subroutine. By establishing the correspondence between data structures in this 
way, the interface effectively sets-up the system of equations to be solved. 


Once the interface has established a correspondence between the numerics’ 
and the physics’ data structures, it is able to perform its third function: providing 
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communications between the numerics and those of the physics. In order to solve 
the system (1.1)1.4), the integrator must be able to evaluate the functions F, 
G, H, and J for different values of U, V, and W. The interface facilitates 
this evaluation using the correspondence which it has established between the 
data structures of the numerics and the physics. To evaluate Ff’, G, H, and J, 
the interface calls the appropriate physical subroutines, passing the appropriate 
values from U, V, and W, and supplying parameter values where necessary. 


1.2.4.2 The Structure of the Interface 


Having described the function of the interface module, its structure is more 
easily explained. As shown in the center section of figure 1, the interface module 
is divided into a control/system set-up section and a communications section. 
As the names suggest, the first of these sections performs the control and system 
set-up functions described above, while the second performs the communications 
function described above. 


The control/system set-up section establishes the system of equations and 
passes this information down along a data pathway (represented in figure 1 
by the dashed arrow) to the communications section. When the program is 
solving the system of equations, control and data flow down along the paths 
represented in figure 1 by the solid arrows, and return in the opposite direc- 
tion. The control/system set-up section calls the integrator to solve the system. 
The integrator, in turn, places a call to the communications section whenever 
it wants to evaluate F', G, H, and J. The communications section, using infor- 
mation which it has received from the control/system set-up section, calls the 
appropriate subroutines with the appropriate variable and parameter values. 
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Chapter 2 


The Physics Module 


The physics module comprises the model’s “knowledge” of the physical pro- 
cesses occurring in a fire. It consists of a set of subroutines which compute the 
functions associated with the algebraic and differential equations describing these 
processes. These subroutines are called physical subroutines and the equations 
which they embody are called physical equations. 


This section describes the forms of physical equations and physicai sub- 
routines and shows how a physical subroutine implements a physical equation. 
It should therefore be of interest to the physicist or programmer who wishes to 
add new equations or subroutines to the model. The section does not include 
a step-by-step procedure for introducing new equations or subroutines into the 
model because understanding such a procedure requires some knowledge of the 
structure of the interface as well as that of the physics. Such a procedure is 
outlined in Appendix A. However, the user may need to refer to this section for 
clarification when attempting to introduce a new equation. Note that this sec- 
tion does not describe the particular physical equations or subroutines of CFC VI 
except by way of example. 


§2.1 The Physical Equations 


2.1.1 A Single Equation in the System of Equations 


The physical equations are a set of algebraic and differential equations which 
embody the physicist’s description of the several processes occurring in the fire. 
Each physical equation states some relationship among the variables in the 
model. 
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Given one physical equation for each currently relevant variable, one hopes to 
completely determine the values of the variables as a function of time. This set 
of equations is called the system of equations. The discussion of the system of 
equations is saved for the sections on the numerics and interface. (The numerics 
section specifies the form of the system of equations and discusses the properties 
which it must have to be solved efficiently. The interface section describes how 
the collection of equations from the physics is bound into such a system.) 


For the moment, simply note that the collection of all the physical equations 
does not, in itself, form a system of equations. First, there may be more physical 
equations than there are variables. For example, if a piece of initially cold P.U. 
foam is ignited with a match, one calculates the spread of the fire by 


ee ele cereale (2.1) 
dt oTs 20T+  307T% ; 


where Fy is the radius of the fire, A is the fire spread parameter, ¢ is the net 
flux impinging on the surface, o is the Stefan-Boltzman constant, and 7’ is the 
flame temperature [4]. But the fire radius of a burning pool is given by 


Ry = Frn(1 a c(t) 2), (2.2) 


where R,, is the maximum burning radius, ¢t is time, and tj is the time of 
ignition [4]. Only one of these physical equations is applicable to a particular 
fire at a particular time in the run, hence only one enters into the system of 
equations for that object and time. There are also cases where only one equation 
exists to describe a particular process, but does not always apply. For example, 
the height of the flames over a burning object is given by the equation 


Ryn 


= ——_—_—_—_, (2.3) 
X|mz| tan Yo 


Hy 


where H, is the height of the flames, R yf is the radius of the fire, 7m, is the burning 
rate of the fuel, Yo is the semi-apex angle of the cone modelling the flames of the 
object in free-burn, x is the fraction of combustion energy released in free-burn, 
and my is the mass loss rate of the object (see subroutine HEIGHT). But this 
equation applies only for burning objects. If the object is not burning, ry is 
0 and the above equation is not applicable. One must then fix the FORTRAN 
variable representing Hy at 0 or, alternatively, leave it undefined. 
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Thus, the collection of physical equations must be properly interpreted before 
it may be thought of as a system of equations. One must find some way to specify 
and choose the appropriate equations to form a system at any given time. The 
following discussion looks at some of the ways in which the system might be 
specified and constructed. 


Some of the difficulties illustrated above could be handled within the physics 
module. Thus in the second example, one could define a logical parameter 
BURN for each object which would be set . TRUE. if an object were burning, and 
.FALSE. otherwise. Then one could write 


__ JReriy/(|\myz|x tan Yo), if BURN is . TRUE.; 
‘iia Ne otherwise. (2st) 


Equation (2.4) is like (2.8) except that it depends on the additional parameter 
BURN and is applicable to all objects, burning or not. The one advantage to this 
approach is that all of the knowledge of the physics is kept within the physical 
subroutines. But the approach runs into difficulty when one tries to apply it to 
a situation like that given in the first example. The numerics could not handle 
a relationship for a variable which is sometimes given by a differential equation 
and sometimes by an algebraic equation. Moreover, in cases like that illustrated 
in the second example (a fairly frequent occurrence) one does not really want 
to evaluate the function Hy at all for an object that is not burning, as that 
would be a waste. One can achieve time savings by fixing Hy at 0 and no longer 
thinking of it as a variable so long as the object is not burning. If this procedure 
is generalized to other variables which are fixed for a portion of the run, the 
time savings can be considerable. Finally, in lumping all the equations for a 
given variable into a single equation with several different cases, one makes each 
variable appear to be dependent on more input variables than it actually is at 
any given time. This can lead to inefficiencies in variable elimination (see section 
4.2.3.1). 


For the reasons given, the problem of multiple equations for the same variable 
is generally handled in the interface rather than in the physics. The interface 
gets information on the state of the fire and uses it to choose the appropriate 
set of equations for the system at a given time. One therefore says that the 
interface has the job of binding the individual physical equations into a system 
of equations. 


The user should note that nothing prevents him or her from writing equations 
(like (2.4) above) which handle several different cases for the same equation within 
the physics. At worst, such an equation will introduce some small inefficiency 
into the program. (S)he should be aware, however, that the apparatus for 
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removing these inefficiencies exists within the interface if (s)he should wish to 
take advantage of it. Thus, a physical subroutine may be applicable to all fires 
or may be restricted to a certain subset of special cases. 


2.1.2 The Forms of a Physical Equation 


A physical equation must have one of three forms if it is to be solved by the 
CFC VI numerical package: 


case (a) 0 
case (b) du,/dt » = $(p,v, t) (2.5) 
case (c) Us 


where p = (p;,..., Dm) are the physical parameters, vy = (v1,...,Un,) are the 
physical variables, ¢ is the time since the start of the run, and ¢:R”t"+!4R 
is the function associated with the physical equation. (Here, IR stands for the 
real numbers, and the notation says that “¢ is a function that maps m-+n-+1 
_ tuples of real numbers into the real numbers.” ) Note that physical variables here 

refers to all those inputs to physical subroutines for which the numerical package 
seeks solutions. Physical parameters refers to all the other inputs to physical 
subroutines except time, ¢. Thus, “physical parameters” need not refer only to 
constants of nature like 7 or the gravitational constant. Dimensions of rooms 
and heights of objects are also physical parameters. Also note that although 
t appears in equations (2.5), in principle actual physical equations should be 
independent of t. In practice, t is used in certain equations that are only crude 
approximations to reality. For example, the fire radius of a pool fire is allowed 
to build up slowly as a function of the time since its ignition in order to ensure 
numerical stability. — 


In case (a) of equation (2.5), the equation is algebraic and written in the most 
general possible form. One calls such an equation a root-finding equation because 
the solution of the system must be a root of the equation. Case (c) is another form 
of algebraic equation. One calls an equation in this form a fixed-point equation 
because the function value it produces at a solution of the system is the same as 
the value of one of the components of that sclution. One can evidently convert 
(c) to the form of (a) by subtracting v, from both sides of the equation. However, 
it is not always possible to solve an equation in the form of (a) analytically for one 
of the variables to get the form of (c). In general, one can solve a system of fixed- 
point equations more easily than a system of root-finding equations (provided 
that each variable in the system appears on the left hand side of exactly one 
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of the equations) since fixed-point equations yield to variable elimination and 
may be solved by a successive-substitution technique (see sections 4.2.3.1 and 
3.1.3, respectively). Therefore, one should try to write algebraic equations in 
fixed-point form whenever practical. Currently, the only algebraic equations 
written in root-finding form are the mass-conservation equations (see subroutine 
TMNET). 


Case (b) of equation (2.5) is the only allowable form for a differential equation 
to be solved by the current numerics. Any ordinary differential equation in time 
can, however, be expressed in that form as follows: Any first-order differential 
equation can be written in the form 


| dv;, dv; du; 
0= OC; Mya» Var) Ups ey —F, «) 5) (2.6) 
where {v;,,...,U,,} are the variables. To express (2.6) in one of the forms of 


(2.5), introduce new variables w,,...,w, and replace (2.6) by the set of g + 1 
simultaneous equations | 


O = (DP; U5.) --+, Uys Wi, +++) Wy; t) (2.7) 
du,;, /dt= wy; 
dv;,/dt= we 

(2.8) 
du, /dt= w, 


Equations (2.7) and (2.8) together have the same solutions as (2.6). Moreover, 
(2.7) is a root-finding equation like case (a) of (2.5) and each equation of (2.8) 
is a differential equation like case (b) of (2.5). A similar procedure converts an 
n‘4.order differential equation into the forms of (2.5). 


As an example, consider the differential equation 


dm 

Hy = Ry src (2.9) 
l-a4|x tan Yo 

which gives the height H; of the flame over a burning object. Here wp is the semi- 

apex angle of the flames in free burn, x is the fraction of energy of combustion 

released in free burn, Ay is the fire radius, ms is the total mass of fuel burned, 

and m; is the total mass loss of the object since ignition. Expressed in the form 
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of (2.6), this becomes 


amy dmyr 


eT eT D ) (2.10) 


0= $(v0,x,Hy, Ry, 


where ¢ is defined to make (2.9) and (2.10) equivalent. Then define new variables 
my, the burning rate of the fuel, and mm;, the mass pyrolysis rate, and replace 
(2.10) by the three equations 


0= O(W0,x, Hs, Ry, mo, my) 


dmy Hn 
—SEE b 
dt (2.11) 
dm; } 
—— = Vy} 
at f 


In practice, the necessity of rewriting the differential equations so that the 
derivatives appear only on the left hand side of the equation presents only a 
minor inconvenience. 


On the other hand, the numerical package used to solve the algebraic and 
differential equations could be easily generalized to handle differential equations 
in the form of (2.6) directly [3]. This project may be undertaken in the future 
though it does not seem to be worth the effort at present. 


2.1.3 Restrictions on the Physical Equations 


Each physical equation must satisfy a number of continuity conditions if it 
is to be solved efficiently with the present numerical package (or by most other 
numerical packages as well). Most of these are consequences of the restrictions 
on the system of equations given in Appendix B, though not every restriction on 
the system of equations can be phrased in terms of restrictions on the individual 
equations. This section presents some of these conditions. 


Given a physical equation of the form (2.5) fix the parameters at p = c where 
c is a vector of constants. The result has the form 


0 
dus ja = Sel, t) (2.12) 


Us 


where f,(v, t)=¢(c,v,¢) defines f,. Restrictions on the physical equations are 
phrased in terms of restrictions on fy. 


14 2. The Physics Module 


First and most crucial, f, must be continuous in all variables and int in some 
neighborhood of the solution. Failure to heed this restriction usually causes the 
numerics to fail at the discontinuity. Moreover, since a computer can only carry 
out calculations to a finite precision, f, begins to look discontinuous when any 
of its partial derivatives become too large. Therefore, one should restate the 
restriction as 


Oke eu (2.13) 
Ou; 


for every variable v; and also 


Ofc 

M 
at 
both in some sufficiently large neighborhood of the solution. It is understood in 
the above restriction that the partial derivatives are to be computed by finite 
differences so that the restriction makes sense even where f, is not differentiable. 


In practice, it is not very difficult to structure the physical equations so as 
to avoid the above type of “discontinuity.” However, there is one subtlety in 
the above restriction which has caused trouble. Suppose one tries to use the 
(hypothetical) equation 


A >he: 


(1/3)chymr7, otherwise; 


to compute the flux, F’, from the flames of a burning object to the surface of that 
same object, where ho is the height of the object, h,; is the height of the interface 
between the cold and hot layers, ry is the fire radius, c is a constant, and hy 
is the height of the flames over the object. Thus, the variables in the equation 
are F’, h;, rf, and hy, and the equation says that the flux is 0 if the object 
is immersed in the hot, sooty layer and is otherwise proportional to the flame 
volume. One might reason as follows: “As the layer interface descends toward 
the surface of the object, the flames of the object receive less and less oxygen 
and the flame height therefore approaches 0. Hence the flux should not change 
discontinuously as the hot layer covers the object.” Though this reasoning is 
quite correct, the above equation may nevertheless cause the integrator to fail as 
the layer covers the object. Although F' does not change discontinuously along 
the path of the solution, the function which computes F is discontinuous in ) f 
at points arbitrarily close to the solution path. This type of discontinuity can 
cause problems in trying to estimate partial derivatives for the Jacobian matrix 
or if small truncation errors accumulate in the variables. To avoid the above 
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discontinuity, rewrite (2.14) as 


0, if ho Shy; 
Cs tC /3)chynr?, otherwise. (2.15) 
For an example of how another type of discontinuity can be removed, see how 
the function ROUND removes a singularity from the equations in FLOW. 


A second restriction on the physical equations is much more difficult to 
satisfy, but fortunately is less crucial. Namely, in order to guarantee a solu- 
tion, the numerics package requires that the function associated with each physi- 
cal equation have continuous first derivatives. In practice, it is nearly impos- 
sible to remove the slope discontinuities or “kinks” from the physical equations. 
However, the numerics routinely obtains solutions in the neighborhoods of such 
kinks, though sometimes with difficulty. At no time has a discontinuous deriva- 
tive been clearly implicated as the culprit in preventing the numerics from finding 
a solution. But since the numerics cannot guarantee a solution in the neighbor- 
hood of a kink, the user would be well advised to avoid writing them into physical 
subroutines where possible. 


§2.2 The Physical Subroutines 


The physical subroutines are the program’s implementation of the physical 
equations. Each physical subroutine computes one or more of the functions 
associated with the physical equations. It takes input parameters and variables 
and produces output|s] which represent the value(s] of the function(s] associated 
with the given physical equations]. Each physical subroutine must conform to 
a fairly strict coding standard. In general, it must be functional: its outputs 
must depend only on its current inputs and not on any previous inputs. A 
few exceptions to this rule are described at the end of this section. It must be 
modular: it may have its own data structures, but must not assume any particular 
data structure in any other section of the program. All communication with the 
interface must be through the argument lists. This section describes the above 
restrictions in detail and suggests a format for coding physical subroutines. 


2.2.1 Physical Subroutines and Physical Equations 


Each physical subroutine evaluates the function associated with a physical 
equation or the functions associated with several physical equations. Recall that 
each physical equation has one of the forms given in (2.5). There, the function 
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SUBROUTINE NETMAS(CP,ZKAZZ,TELZZ, TMDSM, TMLZZ, TMOSM, TMUSM, 
| TMNET1) 
IMPLICIT REAL*8 (A-H,O-Z2) 

CHECKS FOR CONSERVATION OF MASS IN A GIVEN ROOM. 

EQUATIONS BY H. MITLER. CODED BY J. GAHM, SEPTEMBER 1981. 


INPUT PARAMETERS: 
CP -- SPECIFIC HEAT OF AMBIENT AIR. 
ZKAZZ -- AMBIENT TEMPERATURE OF THE ROOM. 


INPUT VARIABLES! 

TELZZ -- RATE OF INCREASE OF THE ENERGY OF THE HOT LAYER. 

TMDSM -- TOTAL MASS OUTFLOW FROM THE COLD LAYER THROUGH 
ALL THE VENTS IN THE ROOM. 

TMLZZ -- RATE OF INCREASE OF MASS OF THE HOT LAYER. 

TMOSM -- NEGATIVE OF SUM OF MASS PYROLYSIS RATES OF ALL 
OBJECTS IN THE ROOM. 

TMUSM -- TOTAL MASS OUTFLOW FROM THE HOT LAYER THROUGH ALL 
THE VENTS IN THE ROOM. 


OUTPUT VARIABLE: 
TMNETI1 -- NET RATE OF INCREASE OF MASS WHICH HAS NOT BEEN 
ACCOUNTED FOR. NOTE THAT WHEN THE SYSTEM HAS 
CONVERGED, WE SHOULD HAVE TMNET#@. 
TMNET1 = TMLZZ = TELZZ/(CP*ZKAZZ) + TMUSM + TMDSM + TMOSM 


RETURN 
END 


a aanannaaaadaaaaaaaaaaagngaana 


Figure 2. Subroutine NETMAS. 


associated with the equation is ¢. A physical subroutine takes its input from 
among the physical parameters, physical variables, and time. It produces one or 
more outputs, each corresponding to the value of some function like ¢, above, 
evaluated at the input point. Given in figure 2 is a typical example of a physical 
subroutine, NETMAS. It takes two input parameters, ZKAZZ (the temperature, Ty, 
of the cold layer in a room) and CP (the specific heat, cy, of air) and five input 
variables, TELZZ (the rate, H,, of increase of energy of the hot layer), TMDSM (the 
net rate of mass outflow, 7ra, from the cold layer of the room through the vents), 
TMLZZ (the rate of increase, mz, of mass of the hot layer), TMOSM (the negative of 
the net mass pyrolysis rate, my, of all the objects in the room), TMUSM (the net 
rate of mass outflow, ™,, from the hot layer through the vents). It computes 
one output variable TMNET (the mass appearing in the room which has not been 
accounted for). 


Note that when the system has converged, one must have TMNET= 0, ie. mass 
must be conserved in the room. Thus, NETMAS computes the function associated 
with the root-finding equation 


Ey 
0=>=my, — 
J Cpl, 


+ My + Ma + my (2.16) 
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Each physical subroutine must evaluate the function associated with at least 
one physical equation. It may evaluate more than one function as well. First, 
since a building may contain several rooms, vents and objects, a physical sub- 
routine may be called several times with different parameters to compute cor- 
responding function values for different entities; e.g. subroutine HEIGHT is called 
once for each burning object to compute the flame height for that object. Second, 
a single call to a physical subroutine may return several different function values; 
e.g. a single call to FLOW returns the mass flow of hot gases, TMUZZ1, through 
a vent, the mass flow of cold gases, TMDZZi, through the vent, and the energy 
flow of hot gases, TEUZZ1, through the vent. In general, subroutines which 
return more than one function value per call can lead to inefficiencies in variable 
elimination (see section 4.2.3.1) and should therefore be avoided where possible. 


Physical equations which reference the data structures present a major dif- 
ficulty for the structure of the program. Consider, for example, the physical 
equation 


>) thu(V) = thu(R) (2.17) 


allVinR 


which just says that the net mass outflow from the hot layer in room R is 
the sum of the mass outflows through all the vents, V, in room R. A physical 
subroutine written to evaluate the left hand side of this equation must receive 
a list of all the vents in room R with corresponding mass flows as input. The 
overhead for the manipulations necessary to create and pass such data structures 
in FORTRAN vastly outweighs the cost of the simple summation which the 
subroutine performs. Therefore, the current program sacrifices some modularity 
in order to perform summations like the above within the interface instead. (See 
“Data Structures Within a Physical Subroutine,” below, for further discussion 
of this problem.) 


2.2.2 The Structure Of a Physical Subroutine 


This section can serve as a guide for writing a new physical subroutine. 


2.2.2.1 Inputs and Outputs 


All input to and output from a physical subroutine passes through the 
argument list. Inputs may be parameters or variables and are always “read- 
only” arguments—they must never be modified in the course of a subroutine 
call. Define input parameters as inputs whose value can never change during 


18 2, The Physics Module 


the course of a single corrector iteration (see section on numerics). Thus room, 
object and vent geometries and physical constants are input parameters under 
this definition. Also note that time from the start of the run, the burning states 
of objects, and previous time step values are considered to be parameters when 
input to a physical subroutine even though their values will change in the course 
of an entire run. For the value of each is determined only once at the beginning 
of a time step and is not changed for the duration of that time step. Conversely, 
input variables are all those inputs which can change in the course of a single 
corrector iteration. These definitions have been chosen so that the variables are 
exactly those inputs which must be taken into account for variable elimination 
(see section 4.2.3.1). Note that variable elimination does not depend on time, ¢, 
for example. | 


The output arguments are called output variables, though each output is 
actually a function value and not a variable value. Output variables are write- 
only; their content is always ignored on input. 


Input parameters may be of any data type. All input and output variables are 
REAL¥8. Experience has shown that REAL*8 precision is needed to ensure stability 
in the numerics, though REAL*4 precision should suffice when improvements are 
made in the numerics. To maintain the usual FORTRAN naming conventions 
so that real and integer data types need not be explicitly declared, begin each 
physical subroutine with the statement IMPLICIT REAL*8 (A-H,0-Z). Nowhere 
in any of the CFC VI physical subroutines are arrays passed through the argu- 
ment lists. However, array arguments are not forbidden in principle, provided 
that the physical subroutine does not depend on any of the interface data struc- 
tures. (See section 2.2.2.2.) 


Input variables should be given the same name within the subroutine as 
they have in the interface to prevent confusion in variable elimination. Output 
variables for fixed-point equations should be given the same name as the cor- 
responding variable in the interface but suffixed with a “1” to indicate that the 
output is actually a function value and not a variable value. The names of input 
parameters should generally correspond to those of their interface counterparts 
as well. 


2.2.2.2 Data Structures Within a Physical Subroutine 


The physical subroutines of CFC VI were written so as to be “data structure 
free.” No data structures more complicated than scalar variables are passed to 
or from the interface. The rationale behind this design choice is that one wishes 
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to make the physics independent of the interface so that changes in the interface 
data structures do not force changes in the physics. 


However, this choice may have been a mistake for the following reasons: (1) 
the necessity of maintaining correspondence between the names in the interface 
and those in the physical subroutines already violates the independence of the 
physics to a degree; (2) there are certain cases where one is forced to “do physics” 
in the interface to avoid passing data structures to the physics. (See the example 
involving equation (2.17) in section 2.2.1 above.) It may therefore be desirable 
to allow the physical subroutines to use more complicated data structures in 
the future. However, one should insist that the data structures, whatever they 
are, be local to the physical subroutines themselves and that the interface take 
the responsibility for translating between the interface data structures and the 
physics data structures. Unfortunately, the manipulation of data structures is 
so cumbersome in FORTRAN that the two problems mentioned above are likely 
to go unsolved for the near future. 


In order to clarify the relationship between the data structures and the 
physics a higher-level description of the forms which the data structures may 
take is needed. If one wants to write an interface flexible enough to accomadate 
changing data structures, one must first delimit the possibilities for changes. 
Much future work is needed in this area. 


2.2.2.6 Coding Standards for a Physical Subroutine 


This section suggests a format for coding physical subroutines. The user may 
wish to modify this format, but should be careful to include all of the suggested 
documentation for use by future programmers and in variable elimination. 


The format will be illustrated by a specific example. Shown in figure 3 is a 
portion of a typical physical subroutine, FLOW. The elements of the argument list 
are ordered with input parameters first, input variables next and output variables 
last. Most of the existing physical subroutines list the arguments within each 
of these three classes in lexicographic order. The IMPLICIT REAL*8 (A-H,0-Z) 
statement ensures that all floating-point variables implicitly declared under the 
standard FORTRAN naming conventions have the necessary precision. The first 
comment lines of the subroutine give a brief description of its function. The next 
block of comments gives a description of each of the input and output arguments. 
Recall a few points on the arguments of the subroutine: 


a. ‘The inputs are read-only arguments; the subroutine must never modify their 
values. The output variables are write only arguments; the subroutine must 
never use their values as input. 
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SUBROUTINE FLOW(CD,CP,VMAZZ,ZBVZZ, ZHRZZ,ZHTZZ, ZHVZZ, 

1 ZKAZZ1,Z2KAZZ0, ZHLZZI, ZHLZZO, ZKLZZI , ZKLZZ0, ZPRZZI1, 
Z ZPRZZO, TEUZZ1, TMDZZ1, TMUZZ1 ) 

IMPLICIT REAL*8 (A-H,0-2) 


C 
C CALCULATES ENERGY AND MASS FLOW THROUGH A VENT. .. 
week Additional comment lines ommitted. **** 


INPUT PARAMETERS: 

CD -- VENT FLOW COEFFICIENT. 

CP -=- SPECIFIC HEAT OF AIR. 

VMAZZ -- DENSITY OF AMBIENT AIR. 

ZBVZZ -- WIDTH OF THE VENT. 

ZHRZZ -- HEIGHT OF THE ROOM(S) CONTAINING THE VENT (NOTE THAT 
WE ASSUME THAT ANY TWO ROOMS WHICH ARE JOINED BY A 
VENT HAVE THE SAME HEIGHT. ) 

ZHTZZ -- DISTANCE FROM THE CEILING TO THE TRANSOM OF THE VENT. 

ZHVZZ -~ HEIGHT OF THE VENT. 

ZKAZZI -- AMBIENT TEMPERATURE OF THE INNER ROOM. 

ZKAZZO -- AMBIENT TEMPERATURE OF THE OUTER ROOM. 


INPUT VARIABLES: 

ZHLZZI -- DEPTH OF THE LAYER IN THE INNER ROOM (ROOM ON SIDE 1 OF 
THE VENT). 

ZHLZZO -- DEPTH OF THE LAYER IN THE OUTER ROOM (ROOM ON SIDE 2 OF 
THE VENT). 

ZKLZZI -- TEMPERATURE OF THE LAYER IN THE INNER ROOM. 

ZKLZZO -- TEMPERATURE OF THE LAYER IN THE OUTER ROOM. 

ZPRZZI -~- ATMOSPHERIC PRESSURE AT THE FLOOR OF THE INNER ROOM. 

ZPRZZO -- ATMOSPHERIC PRESSURE AT THE FLOOR OF THE OUTER ROOM. 


OUTPUT VARIABLES: 
TEUZZ1 -- ENERGY OUTFLOW FROM THE HOT LAYER OF THE INNER ROOM 
THROUGH THE VENT. 
TMDZZ1 -- MASS OUTFLOW THROUGH THE VENT FROM THE COLD LAYER OF 
THE INNER ROOM INTO THE COLD LAYER OF THE OUTER ROOM. 
TMUZZ1 -- MASS OUTFLOW THROUGH THE VENT FROM THE HOT LAYER OF 
THE INNER ROOM TO THE HOT LAYER OF THE OUTER ROOM. 


NgaAgQANNMNANMNAANNNMNAANANNANANNAAAAAANAANANANANNAANN 


PARAMETER ( G = 9.8DM ) 


“eee The calculations of the subroutine are ommitted here. **** 


TEUC Ze CP Oe TEMP. & TEU 
TMUZZ1 = -TEMP * TMU 
TMOZZ 19 82 =T EMPS Se EMD 
RETURN 

END 


Figure §. Format of a typical physical subroutine. 


b. Any input whose value cannot change in the middle of a time step is to be 
considered an input parameter. Any input whose value can change in the 
middle of a time step is an input variable. 


c. Although the example subroutine returns three outputs at every call, in 
general it is desirable that each physical subroutine return only one output 
for each call. 


Most physical subroutines should be functional: outputs should depend only 
on inputs and not on any previous inputs or any internally stored values. In afew 
special cases, clearly delimited under “Non-Functional Subroutines”, one may 
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- 


IF (output variable .GT. 
1 largest physically reasonable value) THEN 
output variable = largest reasonable value 
ELSE IF (output vartable .LT 
1 smallest physically reasonable value) THEN 
habe’: variable = smallest reasonable value 


Figure 4. The form of a typical limit checking section. 


store internal parameters, but these are exceptional cases. Carelessly written 
non-functional subroutines may (1) violate the assumptions of the numerics and 
therefore lead to incorrect answers and (2) make the code extremely difficult to 
analyze and debug. 


Most of the existing physical subroutines which compute the functions as- 
sociated with fixed-point equations end with a section (see TMP02 below for ex- 
ample) which ensures that the outputs of the subroutine are physically reasonable. 
Such a check usually has the form given in figure 4. 


When writing such a limits routine, however, one should take great care 
to avoid introducing discontinuities into the equations. For example, a flux of 
less than 1000 watts/sq. meter incident on the surface of a heating object is 
physically insignificant in a fire of any reasonable size, and one might therefore 
be tempted to write a “limits check” of the form 


Tee(riur LI. 1000.) flux = 0: 


But this line causes the flux to jump discontinuously between 0 and 1000, and 
may therefore cause the integrator to fail. 


2.2.2.4 Non-Functional Subroutines 


As mentioned above, most of the physical subroutines should be functional; 
that is, their output should depend only on their input and not, for example, on 
any values stored in variables local to the subroutine. In certain circumstances, 
however, it is desirable to perform time integrations within a physical subroutine 
rather than within the numerics, and a subroutine constructed to perform such 
an integration must store some information in local variables. 


To integrate the differential equation 


dv 
= ak eet Ui) (2.18) 
dt 


within a physical subroutine, one might calculate v, at each time step wa 
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Ve = Ue + (F(a ‘ » Vig» to ry h)+ Ses ° .+) Vig, to)) (2.19) 


and thereby convert the differential equation to a fixed-point algebraic equation. 
Here, v7, is the converged value of v, at time tp — h and 2;,,...,0;, are the 
converged values of v;,,..., Ui, Tespectively at time tp — h. The 2; are called 
the previous time step values of the variables and h is called the time step. The 
given integration formula is the trapezoidal rule. 


Notice that the ¥; are “parameters” as defined above since they do not 
change during a time step; they are fixed for the n*" time step once the con- 
verged values at the n — 15* time step have been found. Similarly, the entire 
term f(¥;,,...,Us,,t9 — h) can be considered as a parameter. Thus one could 
evaluate the left hand side of the above equation with a functional physical sub- 
routine requiring the following inputs: parameters Uz, h, t, f(Ui,,..., Vi,, to — A); 
variables v;,,...,U;,- But to use such a subroutine one would have to store 
f(¥i,,...,0i,,to — h) within the interface at the end of every time step though 
it is needed by just one subroutine. The alternative is to write a non-functional 
subroutine and store f(U;,,...,0:,, to — 2) as a local variable. 


One can store such local variables, in any subroutine that requires them, via 
a special call to an ENTRY point at the beginning of the time step. Currently, 
only TMPWi and TMPO2 require such internal storage. Portions of TMPO2 appear 
in figure 5. 


At the beginning of each time step a call for each object is placed to TMPO2T 
to update the discrete grid temperature profile for that object. But during the 
time step all calls go to the TMP02 entry. A call to TMPO2 calculates a value for 
ZKOZZ1 using only the input arguments and values in the arrays ZKO, NOB, AO, 
and BO. A call to TMPO2I calculates new values for the elements of the above 
arrays using only old values from ZKO and values from the TMPO2I argument list. 
Therefore, between calls to TMPO2I, calls to TMPO2 are functional; no internal 
parameters in the subroutine change when the TMPO2 entry is called. 


Any physical subroutine which stores internal parameters should conform to 
a format similar to that of TMPO. In particular, it is crucial that every physical 
subroutine be functional within a given time step; ie. that the outputs depend 
only on the input variables within a given time step once the parameters are 
fixed. 


TMPW1 and TMPO2 perform internal integrations because each of these sub- 
routines solves the heat equation wa a numerically calculated double integral 
which would be difficult to remove from the subroutine. In principle, one could 
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SUBROUTINE TMPO2(DT,EB,KO,2ZHOZZ,ZHRZZ, 
i ZKAZZ,Z0AZZ,FQLOR,FQPOR,FQWOR, 
2 ZHLZZ,ZKLZZ,ZKOZZ,ZKOZZ1 ) 
IMPLICIT REAL*8 (A-H,0-Z) 


CALCULATES THE TEMPERATURE PROFILE WITHIN OBJECT KO USING A DISCRETE 
GRID. THE OBJECT IS HEATED ON THE UPPER SURFACE BY RADIATION AND THE 
HEAT DIFFUSES THROUGH IT BY CONDUCTION. NOTE THAT THIS SUBROUITNE 
DIFFERS FROM THE OTHER PHYSICAL SUBROUITNES IN THAT THE TIME AND 
SPACE INTEGRATIONS NECESSARY TO SOLVE THE HEAT EQUATION ARE PERFORMED 
WITHIN THE SUBROUTINE ITSELF. TO FACILITATE THIS CALCULATION THE 
SUBROUTINE MUST BE INITIALIZED AT THE BEGINNING OF EACH TIME STEP 

VIA A CALL TO THE ENTRY TMPO2I1. 


AAANANAANANANN 


**** Some comment lines and parameter statements ommitted here **** 


THE FOLLOWING ARRAYS ARE NEEDED TO PERFORM THE INTEGRATION WITHIN 
THE PHYSICAL SUBROUTINE. 


DIMENSION FIRST(MO1),NOB(MO1),A0(MO1),BO(MO1),ZKO(MXLYRS,MO1 ) 
THE FOLLOWING LINE ENSURES THAT THE DATA IN THESE ARRAYS WILL NOT 


BE LOST UPON EXIT FROM THE SUBROUTINE. 
SAVE FIRST,NOB,AO,BO,ZKO 


QAngna aan 


**** Calculations of the main body of the subroutine ommitted. **** 


Cc 
C LIMITS: 
IF (€ 2ZKOZZ1 .GT. 2828.08 ) THEN 
ZKOZZ1 = 2888.09 
ELSE IF ( ZKOZZ1 .LT. ZKAZZ ) THEN 
ZKOZZ1 = ZKAZZ 
END IF 


RETURN 


c 
C 

ENTRY TMPO2ZI(KO,2G0Z2Z,DTMAX,ZNOZZ,2J0ZZ,2ZKAZZ,DT,ZKO2ZZ) 
c 


C CALL THIS ENTRY ONCE FOR EACH OBJECT TO INITIALIZE AT THE 


C BEGINNING OF EACH TIME-STEP. 
C 


¥*e**e Calculations and comment lines ommitted here. **** 


RETURN 
END 


Figure 5. Example of a non-functional subroutine. 


use this entry point implementation of the trapezoidal rule to convert all of 
the differential equations in the system to fixed-point equations. CFC V, the 
predecessor of CFC VI, actually does use this approach in order to simplify 
the numerical problem. (One needs more sophisticated techniques to solve a 
system of algebraic and differential equations than to solve a set of fixed-point, 
algebraic equations.) However, integrations performed within the physical sub- 
routines may cause errors to accumulate in the results leading to instabilities 
and incorrect answers. As CFC VI provides the numerical machinery necessary 
to solve a mixed algebraic/differential system, the user should not attempt to 
perform integrations within the physics except where necessary. 
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§2.3 Thoughts for Future Development 


The main weakness of the current organization of the physics is that the 
physicist must spread his or her knowledge over too many sections of the program 
in order to incorporate a new physical equation. Ideally, all of the physical 
knowledge incorporated into the program should be conceptually localized in a 
few easily accessible locations. At the moment, only the physical subroutines 
are localized in this way. In order to incorporate a new physical equation 
into the existing program, the programmer must, in addition to writing a new 
physical subroutine: (1) specify the circumstances under which the new equation 
applies and incorporate that knowledge into the interface; (2) possibly perform 
summations within the interface; (3) update the variable elimination scheme; 
and (4) specify how the subroutine interacts with the interface data structure. 
(For an example of how complicated this last problem can be, see the section of 
EVALFP which calculates the vent flows.) When these problems are solved, the 
guidelines for “Introducing a New Physical Subroutine” (Appendix A) can be 
made into an algorithm. 


Chapter 3 


The Numerics Module 


By appropriate manipulation, the physical equations which describe a par- 
ticular fire during a particular portion of its progress for purposes of the model, 
may be written in the form 


0 = F(U,V,W,?) (3.1) 

dV /dt = G(U,V, W,t) (3.2) 
W = H(U,V,W,t) (3.3) 

Y = 1(U,Y,W,t) (3.4) 


where dimn F = dimnU, dimnG = dimnV, dimn H = dimn W, and dimn Y = 
dimn/. The numerical package solves such a system of equations, giving values 
for U, V; W, and Y as a function of time, t. 


The numerical package used in CFC VI is an extension of the algebraic and 
differential equation solver designed by C. W. Gear [3]. It is a variable order, 
variable time-step integrator which chooses the order of the method and the step 
size based on estimates of the local truncation error. 


_ For the most part, this document treats Gear’s integrator as a “black box.” 
The above reference gives a theoretical description of Gear’s algorithm. This 
document merely describes those aspects of the integrator peculiar to the CFC VI 
implementation. 


§3.1 Subroutine ALDIFF 


ALDIFF is CFC VI’s implementation of the integrator. 
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3.1.1 Starting the Integration 


To begin the integration of a new system of equations, program MARK6 places 
a call to entry ADINIT of subroutine ALDIFF. This call sets up certain local 
variables in ALDIFF which it will need in the course of a normal integration cycle 
and resets the order of the method to one. The call to this entry must pass the 
following inputs: 


x the vector (U,V,W,Y) of converged values for the variables 
at the time stored in T; i.e. a solution to the new system of 
equations; 


TIME the time (since the start of the run) at which the new system 
of equations becomes applicable; 


HSTORE the size of the time step to which the estimates of the derivatives 
will be scaled on output (see the expression for X given below); 


EPS the maximum allowable one-step truncation error for the in- 
tegration; 

FNUM _— the number of root finding equations in the new system; i.e. the 
dimension of F’ and U; 


GNUM the number of differential equations in the new system; i.e. the 
dimension. of G and V; 


HNUM the number of fixed-point equations in the new system; i.e. the 
dimension of H and W; 


INUM the number of eliminated variable equations in the new system; 
i.e. the dimension of J and Y. 


On output, ADINIT returns: 


X the values of the variables passed as input, followed by an 
estimate of the first time derivatives of these variables scaled 
by a factor of HSTORE; 


ORDER 1; 


and ALDIFF is prepared to begin normal integration of the new system of equa- 
tions. 


3.1.2 A Normal Integration Cycle 


To integrate the system (3.1)-(3.4), program MARK6 places a call to the 
primary entry, ALDIFF, of the subroutine by that name. A call to ALDIFF must 
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pass the following inputs: 


ORDER 
TIME 


H 
HSTORE 
X 


EPS 
MAXVRS 
MAXODR 
ERROR 


the order of the method to be used in the upcoming time step; 


the time at which the variables in X converged; i.e. the end 
time of the previous time step; 


the largest time step which may be taken during this step; 
the time step to which the derivatives in X are scaled; 


the vector X, given by 
x= (Uttn) V(tn), W(tn), Y (tn); 


heU (ta), PeV (te), ReW (tn), BeY (ta), «- 5 


nk nk hk Bee 
8 rk) _§ 17(k) yi) _s y(k) 
k! U (tn), kt V (tn), k! (tn); kt! Y (t)) 


where & is the order of the method stored in ORDER, h, is 
HSTORE, ¢,, is TIME, and a superscript of “(i)” means “i*# time 
derivative;” 


the magnitude of a significant value for each variable (corres- 
ponding to the storage order in X); 


the smallest allowable time step; 

the maximum one-step truncation error; 
the declared row dimension of JAC; 

the largest allowable order for integration; 


the values output to this array in the previous call to ALDIFF; 


On output, ALDIFF returns: 


X 


TIME 


the latest converged values for the variables and their scaled 
derivatives (i.e. the value of the vector X at time tn+1), or 
the values that were passed as input if the integrator failed to 
complete a step; 


the time at which the integrator converged (ie. the value of 
tn4+1), or the input value if the integrator did not converge; 


a suggestion for the size of the next time step; 
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RESULT a non-zero integer which has a factor of 2’ if the local trunca- 
tion error was too large N times, a factor of 3” if the corrector 
failed to converge M times, a factor of 5 if the Jacobian matrix 
was singular and a factor of —1 if the integration failed; 


ORDER the order of integration for the next time step; 


OLDX the values which were contained in X on input. 


3.1.3 Notes on the Internal Workings of ALDIFF 


During a predictor-corrector cycle, ALDIFF reduces the integration problem 
to that of solving a system of non-linear, algebraic equations. Solving this system 
of equations is called correction and is handled by subroutine CORECT, which is 
designed to find solutions using a multi-variate Newton’s method. 


However, CORECT requires an estimate of the Jacobian matrix, stored in JAC 
in U/L decomposed form, as input and, by controlling the estimation of the 
Jacobian, ALDIFF controls the method used for correction. First, ALDIFF does 
not estimate the Jacobian every time step, but only when the corrector fails to 
converge or when 30 corrections have elapsed since the last Jacobian estimation. 
This saves computation time since decomposing the full Jacobian Matrix is an 
expensive procedure. 


Even when ALDIFF estimates the Jacobian (using subroutine EVJAC), it does 
not estimate the full Jacobian. First, the structure of the equations implies 
that the Jacobian columns for the eliminated variables are known a prion. 
(See section 4.2.3.1 on variable elimination.) They are 0 except for 1’s on the 
main diagonal. Therefore, ALDIFF need not apply the U/L decompostion to 
the full Jacobian, but only to the upper left hand corner square matrix whose 
dimension is given by the number of non-eliminated variables in the system; i.e. 
dimn U + dimn V + dimn W. When CORECT solves the full linear system needed 
for Newton’s method, it first solves for the non-eliminated variables using the 
U/L decomposed portion of the Jacobian, then uses back-substitution with the 
remainder of the Jacobian to find values for the eliminated variables. Here is 
where the benefits of variable elimination are realized. 


Second, ALDIFF often does not even evaluate the Jacobian for all of the non- 
eliminated variables. When the off-diagonal entries of the Jacobian are small, 
one can often obtain convergence in the corrector by assuming that they are 0. 
In this case, the Jacobian is estimated by assuming it to be the identity matrix 
except in those columns corresponding to root-finding variables. For root-finding 
variables, the assumption that the off-diagonal entries are small compared to the 
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diagonal entries is seldom valid. Therefore, these columns of the Jacobian are 
estimated in the usual manner. In this case, only the upper left-hand corner 
square matrix of dimension dimn F’ need be subject to U/L decompostion. All 
of the remaining variables can be found by back-substitution. Consequently, 
this technique is by far the least expensive correction technique. Note that 
when this simplified Jacobian is used, the resulting correction technique is a 
successive-substitution-like technique. 


The current version of ALDIFF switches back and forth between the Newton- 
like technique and the successive-substitution-like technique. It uses the faster 
successive-substitution technique until the corrector fails to converge at twice 
the minimum time step and then switches over to the more robust method until 
the faster method begins to work again. 


The current program’s implementation of the switching between the two 
techniques could probably be improved upon. Since the faster method is less 
stable, the integrator tends to take more time steps when it is in use, cancelling 
most of the savings from the smaller Jacobian decompostions. Moreover, it is not 
clear that the successive substitution technique will work for systems of equations 
which include other root-finding equations; the upper left corner sub-matrix may 
become singular. 


§3.2 Restrictions on the System of Equations 


The integrator (and most, if not all, other integrators as well) cannot guar- 
antee a solution to a system of equations unless the system satisfies certain 
conditions. These have been outlined by J. D. Ramsdell in a note included as 
Appendix B. 
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Chapter 4 


The Interface Module 


The interface module provides communication between the numerics and the 
physics and controls the integrator and the system of equations. Figure 1, in 
chapter 1, is a block diagram of the relations among the two main sections of 
the interface, the numerics, and physics modules. During a typical run, the 
program spends most of its time passing information and control along the path 
of the solid arrows. Suppose that at time to the program is solving the system 
of equations 


0 = Fo(U,V, W,t) 
dV /dt = Go(U, V, W, t) 
W = H(U,V,W,t) (41) 
Y =1,(U,V,W,t). 


(The integrator is designed to solve a system of equations in this form. See 
chapter 3.) Having obtained a solution (U(to), V(to), W(to), Y(to)) of the system 
at time to, the control section of the interface calls the integrator to find a 
solution to (4.1) at time tp + h. In order to find this solution, the integrator 
must be able to evaluate the functions Fy, Go, Ho, and Jo. It evaluates these 
functions vtaa call to the communications section of the interface. 


The communications section of the interface receives values for the com- 
ponents of U, V, and W, and for ¢ as input and provides the function values as 
output, obtained via calls to the physical subroutines. In order to perform this 
task, the communications section of the interface must know: 


a. which physical variable corresponds to each component of U, V, and W; 


b. which physical subroutine to call in order to evaluate each of the components 
of Fo, Go, Ho, and Io; 
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c. the order in which the physical subroutines should be called so as to max- 
imize the dimension of J) and minimize the dimension of Ho (see section 
4.2.3.1); 


d. which parameter values to pass to each physical subroutine. 


A specification of these four sets of information is equivalent to a specification 
of the system of equations to be solved. Therefore, as long as the program tries 
to solve any particular system of equations, all this information is fixed. 


If the program were specifically designed to solve one particular system of 
equations, then the above description of the basic information and control flow 
would be nearly complete. The four sets of information listed above could be 
coded directly into the communications section of the interface and fixed for all 
time. But, in fact, the system of equations to be solved changes frequently. In 
particular: 


a. Every run starts with a different set of initial conditions. Therefore the 
parameters passed to the physical subroutines change from run to run. 


b. At every time step, a few of the parameters passed to the physical sub- 
routines (e.g. previous time step values) change slightly. 


c. When a process of the fire is activated or deactivated the set of relevant 
equations and variables may change drastically. For example, when an object 
ignites, its fire radius, flame height, pyrolysis rate, burning rate, etc. become 
relevant variables. 


d. In the future, the program may revise the calling order of the physical 
subroutines whenever the set of equations changes so as to optimize the 
variable elimination scheme. 


Whenever any such change occurs in the system of equations, the control section 
of the interface passes new information to the communications section so that 
the latter knows how to call the physical subroutines to evaluate the new system 
of equations. 


To understand how the control section of the interface sets the communica- 
tions section to evaluate a particular system of equations, consider the following 
analogy. Referring to figure 1, think of the communications section of the in- 
terface as a “black-box” machine with one input pipe (from the integrator), one 
output pipe (back to the integrator), several external mini-boxes (the physical 
subroutines) attached via input and output pipes, and a bunch of knobs and 
switches forming a control panel on top. During a typical integration cycle, 
the communications box receives a message through its input pipe which con- 
tains a set of values for the variables in the current system of equations. The 
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communications box decodes this message, activates each of the mini-boxes, 
representing physical subroutines, in the appropriate order, then encodes the 
results and returns them through the output pipe. The settings on the control 
panel determine the procedure by which the communications box decodes the in- 
put, activates the physical subroutine boxes and encodes the output. Therefore, 
the settings on the control panel determine the system of equations which the 
communications box evaluates. Whenever the control section of the interface 
decides that the system of equations must be modified, it effectively adjusts the 
settings on the control panel of the interface communications box. 


In reality, the five common blocks /VAR/, /BLDNG/, /VARYNG/, /SYSPTR/, 
and /PCKPTR/ perform the function of this metaphoric control panel. The 
values stored in these common blocks correspond to settings of the control panel 
switches and dials. Whenever the control section of the interface determines 
that a new system of equations applies, it modifies these values accordingly. 
Common /VAR/ holds values for any physical variable which is not actually 
varying at that point in the run. Common /BLDNG/ holds a value for each of the 
physical parameters. Common /VARYNG/ tells the communications section which 
physical subroutine to call for each currently relevant variable. Finally, commons 
/SYSPTR/ and /PCKPTR/ tell the communications routines how to interpret the 
vector of variable values passed from the numerics and how to store the output 
information before it is returned to the numerics. 


Thus, whenever the control section of the interface resets the system of equa- 
tions, information passes down through the dashed arrow in figure 1, rather than 
along the usual pathway of solid lines. This alternate information pathway is like 
a hand that reaches down to reset the controls of the hypothetical communica- 
tions box. This is the only place where common blocks are used to transmit 
information between different sections of the program. 


§4.1 The Control Section of the Interface 
Figure 6 shows the structure of the control section of the interface. The 
following discussion describes each of the modules of this section in detail. 


4.1.1 Program MARK6 


Program MARK6 is the main controlling module of the fire code. Figure 7 
shows the structure of MARK6. 


Familiarity with three basic concepts is a prerequisite to understanding the 
structure of this module. The first of these concepts is that of an interruption. 
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open files 
initialize variables 
get input 
while haven't yet reached the stopping time do 
if time for a discontinuity then 
call NWSTAT to get new system of equations 
call ADINIT to reset the integrator to first order 
set TDISC beyond end of run so disc. won't be repeated 
end if 
if time for an output then 
call output package 
note time for next output 
end if 
if time for any other type of interruption then 
act on that interruption 
end if 
set tint to time of next planned interruption 
while time < = TINT do 
call NEWSTP to adjust system of equations 
call SCALE to get values in XMAX 
save enough information to recover in case the 
integration steps beyond a discontinuity 
. get size for the time step 
call ALDIFF to take a step 
if step failed then 
if using a first order method, bomb out 
otherwise try resetting integration to first order 
else (step succeeded) 
call UNPINT to see if integrator stepped beyond disc. 
if so get time of dise. in TDISC 
if disc was detected, set TINT=TDISC 
if disc was detected, recopy all previous time step 
values to back up to beginning of step 
end if 
end while 
end while 
write out final values and stop 


Figure 7. Structure of program MARK6. 


Ordinarily, the program integrates the system of equations, getting converged 
values for the variables at times to, t; == to + ho, to = t; + 4,..., ln = 
tn—1 +h,z—1, but taking no special action at the end of each time step. Under 
these circumstances, the numerical package is free to choose whatever values it 
finds convenient for the h,;’s. It does not matter if the program obtains converged 
values at ¢ = 14.5 and ¢ = 16.5 or at ¢ = 11.8 and ¢ = 16.7. Occasionally, 
however, one must force the integrator to converge at a particular time, ¢;. For 
example, if the user has requested output at time ¢; = 20.0 then the integrator 
must converge at 20.0 seconds so that converged variable values may be output. 
Similarly, if an object reaches its ignition temperature at time ¢;, = 310.56 then 
the integrator should converge at t;, so that the system of equations may be 
properly modified. (The system of equations may not be modified in the middle 
of a time step.) Any event which requires that the program converge at the time 
of its occurrence is an interruption. 


Second, an interruption may be planned or unplanned. It is planned if the 
program can foresee the interruption before the integration reaches the time for 
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Figure 8 Planned and unplanned interruptions and discontinuities. 


its occurrence. For example, if the user requests output to the terminal every 
20.0 seconds of fire time, a planned interruption occurs at t = 0.0, 20.0, 40.0,... 
etc. If, on the other hand, the program cannot foresee the interruption and 
must go past the time of its occurrence in order to detect it, the interruption 
is unplanned. For example, suppose a heating object is set to ignite when its 
surface temperature reaches 940.0 degrees Kelvin. One can examine the surface 
temperature of the object at the end of each time step to see if it has reached the 
critical value. But one cannot predict in advance the exact time at which the 
object will ignite. When the program detects that it has stepped beyond the time 
for an unplanned interruption, it estimates the time at which the interruption 
occurred, backs up to the beginning of the just completed time step, and plans an 
interruption at the estimated time of the event. Thus, unplanned interruptions 
become planned interruptions once they are detected. 


Third, there is a special type of interruption known as a discontinuity. A 
discontinuity is an event which causes a discontinuous change in the values of 
the variables and/or a change in the set of equations to be solved. For example, 
when an object ignites its fire radius jumps discontinuously from 0.0 to its initial 
burning radius (a preset value) and several new equations enter the system. 


Figure 8 illustrates the relationship among the three concepts just introduced. 
Referring to the numbers in the various regions of the figure, an example of 
each type of interruption follows: 


Write out a table of variable values when a temperature sensor reaches a 
certain temperature (not implemented in current program). 


Print out a table of variable values at time t = 210.0. 


Ignite an object when its surface reaches a certain temperature. 


36 4, The Interface Module 


Ignite an object at time 215.0 with a match (not implemented in current 
program). 


Referring back to figure 7, MARK6 consists of an initialization section, two 
nested WHILE control structures (simulated using ordinary FORTRAN control 
statements), and a finishing section. The initialization section opens the output 
file, calls INIT to read the input data file, and initializes local variables. The outer 
WHILE loop is repeated until the integrator reaches the finishing time. It checks 
to see if the integration has reached the time for any planned interruptions and 
takes the appropriate action, locates the next planned interruption and enters 
the inner WHILE loop. Note that when the interruption is a discontinuity, a call 
is placed to ADINIT to restart the integrator at first-order. The inner WHILE 
loop integrates until it encounters an interruption. This may be the planned 
interruption which has just been located in the outer loop or an earlier unplanned 
interruption discovered in the course of the integration. To take a time step, the 
inner WHILE loop calls NEWSTP to set up the system of equations for the new 
time step, calls SCALE to set up XMAX for ALDIFF, saves enough information to 
back up the integrator if it steps beyond an unplanned interruption, selects the 
step size, and calls the integrator. If the integrator fails to converge, MARK6 tries 
to restart it at first order. If this fails, it terminates the run. If the integrator 
successfully completes a time step, MARK6 calls UNPINT to see if any unplanned 
interruptions occurred during the time step. If so, TINT is reset to the time of the 
interruption and the integrator backs up to the beginning of the just completed 
time step. The finishing section merely writes out the final values and stops the 
program. 


4.1.2 Subroutine UNPINT 


Occasionally, in its normal operation, the integrator may accidentally step 
beyond a time at which the program must take some special action. This 
does not happen through error, but because the program cannot detect certain 
interruptions until it has stepped beyond them. For example, the user might 
request output when a temperature sensor reaches 700.0 degrees Kelvin, and 
the program cannot predict in advance the time at which this will occur. Such 
an event is called an unplanned interruption; “unplanned” because the program 
cannot foresee its occurrence and “interruption” because it requires that the 
integrator converge at the time of its occurrence. If an interruption requires a 
change in the system of equations it is called a discontinuity. The description 
of PROGRAM MARK6 above contains a full discussion of planned and unplanned 
interruptions and discontinuities. 
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When the integrator steps beyond an unplanned interruption the program 
must: 


a. detect the interruption and note what type it is; 


b. estimate the time at which the interruption occurred and schedule a planned 
interruption at that time; 


c. back up to the beginning of the previous time step, this time being careful 
not to step beyond the time of the interruption. 


Subroutine UNPINT detects unplanned interruptions and estimates the time 
of their occurrence. 


After each integrator call, MARK6 sends the newly converged variable values 
to UNPINT to check for interruptions. Currently, UNPINT checks for four types 
of interruptions, all of which happen to be a co oan (i.e. they require 
changes in the set of equations): 


1. Burnout of an object occurs when the remaining mass of that object drops 
below the value set by the PARAMETER BRNOTM. 


2. Ignition of object 0 in room R by heating occurs when its surface reaches 
the ignition temperature ZKOIG(0,R). 


3. Ignition of an object by contact with the flames of another object occurs 
when any part of that object lies above the disc of flame on the surface of 
any burning object and below the tip of the cone of these flames. 


4. Spillover into an initially uninvolved room occurs when hot gases from any 
other room begin to flow through the vents into that room. 


When UNPINT detects any one of these discontinuities, it sets TDISC to the 
estimated time of its occurrence, sets CHNGRM to the room number of the room in 
which the event occurred, and sets CHNGOB to the object number of the affected 
object (if relevant). If UNPINT is modified in the future to check for other types 
of discontinuities, it may be required to return other types of information. For 
example, if the physicist wishes to add a provision for opening holes in the ceiling 
of a room when the ceiling reaches a certain temperature, (s)he might require 
UNPINT to return the position and size of the hole, as well as the number of the 
affected room, upon detection of such an event. 


In the current program, the time of occurrence of a discontinuity is estimated 
as follows: 


1. In the event of a burnout, the time at which the mass of the object reached 
the burnout mass is estimated by linear interpolation. 
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2. In the event of an ignition by heating, the time at which the surface tem- 
perature of the object reached its ignition temperature is obtained by linear 
interpolation. 


8. In the event of an ignition by contact with flames, the time of ignition is 
taken to be the time of its detection. : 


4, In the event of a spillover, the time of the discontinuity is taken to be the 
time of the start of the just completed time step. 


At first, estimating the time of occurrence of these discontinuities may seem 
a waste of effort. Since the fire model is generally only a crude approximation to 
reality, why should one care if an object ignites at ¢ = 310.0 or t = 310.5? One 
could just pretend that every discontinuity occurs at the end of the time step in 
which it is detected and thereby save the effort of backing up the integrator to 
converge at the time of each discontinuity. But the machinery for backing up the 
integrator and estimating the time of a discontinuity was not introduced only to 
increase the accuracy of the results. Rather, it is intended to prevent errors and 
instabilities which might arise from a situation like the following: Suppose that 
at time ¢ = 59.5 hot gases begin to flow from room 1 into room 2, which was 
previously uninvolved. The program detects this discontinuity upon converging 
at time ¢ = 60.0 and attempts to restart the integrator without properly backing 
up to t = 58.0, the start time of the just completed time step. The integrator 
might fail for the following two reasons: (1) When the variables pertaining to 
the layer in room 2 are finally drawn into the system at time t = 60.0 the flow 
through the vent has increased to the point where it may cause some of them to 
change discontinuously. If, however, the variables are drawn into the system at 
t = 58.0 seconds, well before the layer spills into room 2, this danger is averted. 
(2) Between ¢t = 59.5 and ¢ = 60, hot gases flowed through the vent connecting 
rooms 1 and 2 but, since the equations pertaining to the hot layer of room 2 had 
not yet been drawn into the system, these gases did not appear in the layer of 
room 2; they were lost. This situation is clearly unphysical and may therefore 
have unpredictable results. 


When more than one interruption occurs within the same time step, UNPINT 
selects the one which is estimated to have occurred first. Once an interruption 
is detected and a time for its occurrence is set in TDISC, no other interruptions 
may occur before TDISC. This ensures that the program does not back up to the 
previous time step, TOLD, step to TDISC and find that the interruption actually 
occurred between TOLD and TDISC, a process that could be repeated indefinitely. 
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4.1.3 Subroutine NWSTAT 


At certain times in the course of a run, the set of algebraic and differential 
equations may change abruptly. At the start of the run, the appropriate equa- 
tions must be drawn into the system and the corresponding variables initialized. 
When an object ignites the equations describing the behavior of its fire enter the 
system, and when it burns out they leave the system. Similarly, when hot gases 
first enter a previously uninvolved room the equations describing the hot layer in 
that room enter the system. Future developments in the physics may require the 
program to handle many more such changes in the set of equations. Subroutine 
NWSTAT updates the system of equations and reinitializes the variables when these 
changes occur. 


Subroutine NWSTAT does not detect the times at which changes in the set of 
equations occur. Program MARK6 and subroutine UNPINT perform this function. 
MARK6 tells NWSTAT when the integrator has reached the time for the equation 
change. NWSTAT receives arguments describing the discontinuity which has been 
detected: DSCTYP gives the type of discontinuity (burnout, ignition by heating, 
ignition by contact with flames, start of the run, or spill-over of hot gases into 
a new room), OB gives the affected object (if relevant), and RM gives the affected 
room (if relevant). 


Given this information, along with the converged solution (contained in X) 
of the old system of equations at the time (stored in TIME) of the discontinuity, 
NWSTAT performs several functions: (1) It returns to MARK6 the number of root- 
finding equations, FNUM; the number of differential equations, GNUM; the number 
of fixed-point equations, HNUM; and the number of eliminated variable equations, 
INUM, in the new system. (2) It returns the solution of the new system at the 
time of the discontinuity. (3) It makes any necessary changes to the physical 
parameters in common /BLDNG/ (e.g. it changes STATECOB,RM) to 5 if object OB 
in room RM ignites at this discontinuity). (4) It sets up the commons /VARYNG/, 
/PCKPTR/, /VAR/, and /SYSPTR/ so that the interface communication section 
evaluates the new system of equations rather than the old. (The auxiliary 
subroutine SETSYS performs this last function. See under that heading for a 
more complete description.) 


NWSTAT performs its task in five steps: 
Call UNPACK to get variables into common /VAR/. 


Reset parameters in common /BLDNG/. 


Call SETSYS to establish the new system to be solved. 


Pe Np 


Get a first guess for the solution of the new system. 


40 4, The Interface Module 


5. Call numerics to get an exact solution for the new system. 


NWSTAT receives the converged variable values in a one-dimensional array, X. 
The call to UNPACK unpacks this vector into common /VAR/ so that each physical 
variable may be referenced by its descriptive name rather than simply as a 
nameless component of a vector. The next block of statements resets the values 
in /BLDNG/ according to the type of discontinuity which has been encountered. 
SETSYS uses these updated values to select the new system of equations. 


The problem that remains is finding a solution to this new system of equations 
based on the solution, now residing in /VAR/, of the old system of equations. The 
integrator requires such a set of initial values for the new problem in order to 
proceed along the solution path of the new set of equations. 


Actually, “initial values” means two different things here. First, one initial 
condition must be specified for every differential equation in the new system. 
This is not just a matter of convenience for the integrator. Every system 
containing n differential equations requires n initial conditions to determine a 
unique solution. Thus, the physicist must be sure that every integrated variable 
is properly initialized at each change in the set of equations. Second, once 
the initial values for the integrated variables are fixed, the program must find 
a set of initial values for the rest of the variables in the system. ‘I’his, in a 
sense, is simply a matter of convenience for the integrator. With values for 
the integrated variables already specified, the initial values for the remaining 
variables are simply the unique solutions of the algebraic equations in the system. 
Once SETSYS sets up the new system of equations, one could theoretically obtain 
the necessary initial values by solving the subsystem of algebraic equations. In 
practice, however, one cannot solve the algebraic equations without a good first 
guess. 


More formally, the problem is this: Suppose a change in the system of 
equations occurs at time tg. Call this discontinuity a. For t < tg the pertinent 
system of equations is 


0 = A;(U,, Vi, t) 
dV; /dt = D,(U,, Vi, t) (4.2) 


with initial conditions 
C1 — Vi (ta), 


where dimn A; = dimn U,; = m, and dimn D; = dimn V, = n,. For t>tg the 
pertinent system of equations is 
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0 = Ao(U2, Vo, t) 


dV2/dt ne D2(U2, Vo, t) (4.3) 
with initial conditions 
Co = Vo(ta), 
where dimn A> = dimnUz = mp2 and dimn Dp = dimnVo = no. Upon 


encountering this discontinuity, the physicist must specify, at the very least, 
a function (peculiar to that discontinuity) 


peas Rm _, pre 


defined so that 
Co = fa(Ui(ta), Vilta)). 


In words, (s)he must provide a rule for computing new initial values for the 
integrated variables from the old values of all the variables. 


- But in order to proceed along the new solution path, the integrator must 
know the complete solution to the new system at time tg. Therefore, NWSTAT 
must find values for U2(tg) as well as Vo(t2). To find these values, one might 
provide a second function 


Da: Rmr1_ppme2 


defined so that 
Uo(ta) = ga(Ui (ta), Vilta)); 


i.e. so that 
0= Aa( ga(Ualta) Vi (ta)), fo(Us(ta), Vi(ta)), ta) 


In fact, this is the best way to obtain the needed values for Ua(tg). But, at 
first glance, it may seem redundant to explicitly construct a function g,. The 
function fa, along with values for U,(tg) and Vi(tg), and the new system of 
equations, completely determines values for the variables of the new system as 
long as that system applies. Why, then, could one not solve the purely algebraic 
system . 


0 = Aog(U2, Co, ta) (4.4) 
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for U2 in order to get the required initial values? In principle, the values for 
U2(tz) could be found in this way without any need for the g, function. In 
practice, however, a system like (4.4) can only be solved by a numerical technique 
which requires a good initial guess. The rub then is that although the physicist 
need not specify a gq if (s)he intends to obtain U2(t,) by solving (4.4), (s)he must 
specify instead a good first guess: a function 


ip (Rmtrm_,pma 


defined so that. 
Ga (Us(ta), Va(ta))—Vatta)]] <«. 


Thus, for every discontinuity, a, the physicist must construct two functions to 
provide initial values for the new set of equations: an fa function to provide 
initial values for the integrated variables and either a g, or a 9, to aid in finding 
initial values for the algebraic variables. These functions are called transition 
functions because they assist in the transition from one system of equations to 
another. 


The physicist reading this section should take care not to skip over the 
above formalism without appreciating its full implications. In practice, it is 
extremely difficult to construct a g-type transition function to provide either an 
exact solution to the new system or an adequate first guess. Moreover, every 
time any modification is made in the physics, the g functions for every type 
of discontinuity must be updated. Thus, the necessity of specifying transition 
functions for every type of discontinuity represents a serious obstacle to the 
flexibility of the program. 


No modification in the interface could do much to ease the difficulty of 
providing transition functions. All numerical techniques for solving systems of 
equations require a good first approximation to the solution. (See Appendix 
B.) However, there are some tricks which the physicist can use to simplify the 
task considerably. First, note that it is generally not too difficult to specify the 
f-type transition functions: At a change in the set of equations, the integrated 
variables usually take on a pre-specified value or do not change at all. For 
example, when an object ignites, its fire radius is set to its initial burning radius, 
a constant, while the mass of the hot layer in the room containing this object, 
another integrated variable, does not change at the instant of ignition. In the 
future, the physicist may wish to introduce instantaneous processes which cause 
more complex changes in the values of the integrated variables. For example, an 
explosion, in a room might cause the composition of the hot layer in that room 
to change dramatically. The f transition function for this discontinuity would 
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calculate the mass of each gaseous species in the layer after the explosion taking 
into account the type, size, and location of the explosion. For the moment, 
however, the f transition functions are always identity or constant functions. 


At present then, constructing the g transition function for a discontinuity is 
the more difficult task. Suppose one tries to construct a g transition function 
for the ignition of an object. The object is ignited with a match at time ¢j,. 
Instantly (according to the model) the disc of fire on the object attains a radius 
of 3.7 cm. Given this new fire, the g transition function must produce at least 
a good estimate for every algebraic variable relevant to the current situation, 
including, for example, the pyrolysis rate of the new fire, and the flux from the 
new flame down to its host object. But the pyrolysis rate depends partly on flux, 
and the flux depends partly on the pyrolysis rate. The interdependencies of the 
other algebraic variables are similarly tangled. In order to construct a g function, 
the physicist must untangle them. As the model becomes more complex, (s)he 
will find this an increasingly difficult task. 


Therefore, the physicist should try to construct the model so that the trans- 
tion functions are as simple as possible. Ideally, they should be the identity 
wherever possible; i.e. the equations should be constructed so that the variable 
values do not change discontinuously even when the system of equations changes. 
Take an example from the existing program. When a real pool of gasoline ig- 
nites, the spread of the flames is so rapid in comparison with the other processes 
of the mode] that one could properly assume that the fire reaches its maximum 
burning radius instantaneously. But deriving a g function for this discontinuity 
would pose a formidable challange. If the pool fire is large, the values of variables 
describing it would be heavily interdependent with the values of other variables 
in the system. Therefore, when a pool fire ignites in CFC VI, the radius of the 
fire starts at 0.0 and is allowed to build up gradually over a period of a few 
seconds of fire time. Then none of the variables change discontinuously and the 
g function is the identity for most variables. The errors introduced by the slow 
start-up are negligible. Most of the transition functions can be simplified by 
similar procedures. 


Using procedures like the above, one could simplify the g transition functions 
to the point where they become exact expressions for the solution of the new 
system of equations. If the solutions of the new system are exactly the same as 
those of the old system at the time of the switch-over, then the identity transtion 
function gives the exact solution of the new system. There is great advantage 
to making g an exact expression in this way. If the exact set of initial values for 
the new system is known a priorz, there is no need to call the numerics to solve 
the algebraic equations by Gauss-Seidel or Newton iteration. Such a system is 
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therefore much more reliable. One need not worry about whether or not the 
numerics will be able to converge at a particular discontinuity. Finally, if one 
uses only an approximate g function, one may find that the function does not 
provide a close enough initial guess for certain problems, even if it works in most 
cases. Therefore, even though the current program uses the numerics to get a set 
of initial values for the algebraic variables given only an approximate g function, 
the proper way to handle the transition functions is to require that the physicist 
give (and therefore be careful to ensure that (s)he is able to give) an exact g 
function for every discontinuity. No call to the numerics would then be needed. 


4.1.4 Subroutine SETSYS 


Whenever the program encounters a discontinuity which requires a change 
in the system of equations to be solved, subroutine SETSYS constructs the new 
system. In particular, SETSYS sets up the common blocks /VARYNG/, /VAR/, 
/PCKPTR/, and /SYSPTR/ so that the communications section of the interface 
evaluates the new system of equations on receiving a call from the numerics. The 
subroutine performs its function in two steps. First, it uses the information in 
/BLDNG/ to determine which variables and equations are relevant to the current 
physical situation and sets up /VARYNG/ and /VAR/ to reflect this information. 
Second, it chooses from among the relevant variables and equations those which 
actually must be sent to the numerics. It stores this information in /PCKPTR/ 
and /SYSPTR/. 


Every fire model to be handled by the CFC VI interface has a maximum set 
of variables constituting the dimensions along which it describes every particular 
fire. (These variables correspond to the elements of section 1 of /VAR/, described 
below.) In relation to a particular fire, each of these possible variables is either 
(1) irrelevant (as in the case of a variable reserved for the description of an entity 
which does not exist in the current run, e.g. the pressure, ZPRZZ(3), in room 
3 during a two room run); (2) relevant but not varying (as in the case of the 
burning rate of a non-burning object); or (3) relevant and varying (as in the case 
of the flame height of a burning object). In the case where a variable is relevant 
and varying, it may vary according to any of several regimes; that is, any one 
of several equations may describe its progress. For example, the current model 
calculates the radius of the disc of fire on an initially cold object ignited with 
a match waa differential equation but uses a fixed-point equation to calculate 
the same variable when the object is a pool of gasoline or is pre-heated to a very 
high temperature before ignition. Therefore, for a given variable, v, at a given 
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point in the calculation, one of the following cases holds: 


(case Oa) v is irrelevant; 

(case Ob) v =c (v is fixed); 

(case 1) equation e; describes the progress of v; 

(case 2) equation é describes the progress of v; (4.5) 
(case n) equation e, describes the progress of v; 


where c is a constant and each e, has one of the forms f(v,...,¥m,t) = 0, 
Gin, tr, t) = G0 /diporh(o;j 0) punters 


The first section of SETSYS looks at the parameters in /BLDNG/ to decide 
which of the above cases applies to each possible variable. The arrays FIROOM, 
LAYR, STATE, FIRTYP, and OBJS and the integers NROOMS, NVENTS, and NWALLS 
are particularly useful in selecting among the cases for each variable. SETSYS 
stores its selection in the common blocks /VARYNG/ and /VAR/. 


SETSYS accesses common /VARYNG/ via the integer array VARY which contains 
one integer element corresponding to each possible variable of the model. For 
a given variable v, SETSYS sets the element of VARY corresponding to v to 0 
if v is irrelevant or fixed, and to 7 if v is varying according to an equation 
é;. In the case where v is fixed at a constant c, SETSYS must also specify the 
value of c. This value is stored in common /VAR/. Common /VAR/ is divided 
into three sections of which only the first is relevant to the current discussion. 
(The other sections of /VAR/ are discussed in the description of EVALFP and 
ERFDE.) The first section of /VAR/ contains one REAL*8 element corresponding 
to each variable of the model in the same way as the elements of /VARYNG/.! 
SETSYS accesses /VAR/ through arrays named for each of the types of variables 


Because of the correspondence between the positions in /VAR/ and the variables in the system, 
there is a tendency to confuse the values in /VAR/ with the values of the variables themselves. 
But /VaR/ is not a storage area for the converged values of the variables. (The one-dimensional 

- array X in PROGRAM MARK6 serves as such a storage area.) It is a working space for the 
convenience of the interface which may or may not contain converged values for the variables 
at a given point in the calculation. Common /VAR/ serves to communicate values for the 
fixed variables from the control section of the interface to the communications section of the 
interface. It should never be used to pass values from the communications section of the 
interface to the control section or to any other part of the program. For example, an output 
package should not simply print out the values it receives in common /VAR/ and expect these 
to represent converged values for the variables. Instead, it should receive the variable values 
in X and unpack them into /VAR/ locally. Failure to heed this advice can make the code very 
difficult to modify, since the programmer will constantly have to keep track of what every 
part of the code does to the contents of /¥VAR/. 
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in the model; e.g. ZRFZZ(0,R) refers to the position in /VAR/ corresponding to 
the radius of the fire on object 0 in room R. Upon determining that variable 
v is fixed at c at a particular point in the calculation, SETSYS assigns c to the 
position of /VAR/ corresponding to v. For example, if object 0 in room R is not 
burning, SETSYS assigns 0.0 to ZRFZZ(0,R). 


Note that it is necessary to specify meaningful values for variables that are 
not varying at a particular point in the calculation because certain physical 
subroutines may take these variables as input. If the value for ZRFZZ(0,R) were 
left. unspecified in the above example, any subroutine which took ZRFZZ as input 
might return garbage. Conversely, note that no output of a physical subroutine 
should depend on the value of any variable which SETSYS considers irrelevant 
((case Oa) of (4.5) above) at a particular point in the calculation. 


Testing each condition in turn, SETSYS selects all of the variables relevant to 
the current problem and sets up an equation corresponding to each. (A variable is 
relevant if and only if (case 0a) of (4.5) does not apply to it.) For a given variable, 
v, if (case Ob) of (4.5) applies, the corresponding equation is v = c. Otherwise, 
if (case i) applies, then the corresponding equation is e;. Thus, the collection of 
all these equations forms a system with as many equations as unknowns. 


By an appropriate ordering of the equations, this system can be written in 
the form 


F"(U,V,W, Z,t) =0 
G"(U,V,W, Z,t) = dV /at 
* 
MCE V, W,, Z, t) a oe (4.6) 
I"(U,V,W,Z,t)=Y° 
C=Z. 


where U = (uj,...,Un) is the vector of currently relevant variables with cor-. 
responding root-finding equations; F” = (f;,...,f,,) is the vector of the func- 
tions associated with these root-finding equations; V = (v1,..., Um) is the vector 
of currently relevant variables with corresponding differential equations; G’ = 
( 935 Pads g,,) is the vector of the functions associated with these differential equa- 
tions; W = (wj,..., Wp) is the vector of currently relevant variables with cor- 
responding fixed-point equations; H’ = (hi, tes hs) is the vector of functions 


associated with these fixed-point equations; Y" = (y},...,y,) is the vector of 
currently relevant variables with corresponding eliminated variable equations; 
I"* = (i;",...,%") is the vector of functions associated with these eliminated 
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variable equations; Z = (z,,...,2,) is the vector of variables that are fixed at 
this point in the calculation; and C = (c,,...,c,) is a constant vector. 


The numerics could solve the system (4.6) directly, but one can obtain greater 
efficiency by removing unneeded variables and equations from the system. First, 
one can remove all the variables in Z from the system to obtain the smaller 
system 


F(U,V,W,t)=0 
G(U,V,W, t) = dV /dt 
H(U,V,W,t)=W 
r(U,V,W,t)=Y" 


where 


F(U,V,W,t) = F"(U,V,W,C,t) 
G(U,V,W, t) = G"(U,V,W,C,t) 
H(U,V,W,t) = H'(U,V,W,C,t) 
I"(U,V,W,t) =I'"(U,V,W,C,t) 


defines F, G, H, and I". Next, one could remove all the eliminated variables 
from the system and instead solve 


F(U,V,W,t)=0 
G(U,V,W, t) = dV /dt 
H(U,V,W,t) = W, 


substituting the solution thus obtained back in J” to find Y". However, it 
happens that some of the eliminated variables are the least stable of the variables 
and it is therefore desirable to send a few eliminated variables to the numerics 
to allow the integrator to check them for convergence and to allow their error 
estimates to control the step size. (The numerics handles eliminated variables 
with very little expense since it does not need to compute their corresponding 
rows in the Jacobian matrix.) Therefore, the system actually solved by the 
numerics is 
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FU, V,W,t)=0 

G(U,V,W,t) = dV /dt 

H(U,V,W,t)=W (4.7) 
I(U, VW; t) =Y, 


where Y is just like Y" but with some of the components removed, and / is 
like J" but with corresponding components removed. In order to choose which 
eliminated variables need not be sent to the numerics, first send all eliminated 
variables to the numerics and then remove those whose error estimates are 
relatively small. 


The second section of SETSYS constructs such a system by setting up pointer 
arrays which define the components of F, G, H, J], U, V, W, and Y. Define a 
vector X =(U,V,W,Y) and a vector of functions B = (F',G,H,J). Then the 
g'“ position of the integer array INSYS in common /SYSPTR/ receives a value 
giving the position of the element of /VAR/ corresponding to the g** component 
of X. Similarly, the g*" position of the integer array PACK in common /PCKPTR/ 
receives a value giving the position of the element of /VAR/ containing the output 
of the g**component of B after a function call. (For a description of how /VAR/ is 
used during a function call, see sections on EVALFP and ERFDE. ) Also, the integers 
FNUM, GNUM, HNUM, and INUM receive values giving the dimensions of F’, G, H, 
and J, respectively. SIZE receives a value giving the total number of equations 
in the system. The auxiliary subroutine PUTPTR aids SETSYS in placing pointers 
into these arrays. 7 


The information in these common blocks and integer variables enables the 
interface to interpret between the interface and numerics data structures so that 
the integrator sees a problem in the abstract form of (4.7) while the physics 
receives calls to concrete subroutines taking physically meaningful inputs. 


4.1.5 Subroutine NEWSTP 


While NWSTAT and SETSYS reset the system of equations whenever it changes 
to include entirely new equations, subroutine NEWSTP makes the minor changes to 
the system of equations which occur at every time step. These include updating 
previous time step values in common /BLDNG/ and calling the initialization 
entries of those physical subroutines which perform internal integrations. Note 
that the changes to the system of equations which NEWSTP makes must not cause 
discontinuous changes in the values of the variables. 
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§4.2 The Interface Communication Section 


In order to solve the system of algebraic and differential equations describing 
the fire, the numerics must be able to evaluate the functions associated with 
that system. In CFC VI, the system of equations has the form given in (4.7) and 
the integrator, ALDIFF, evaluates the functions via a call to the communications 
section of the interface of the form CALL FUNCT(X,Z,T). Subroutine FUNCT 
takes X and T as inputs and produces Z as output, where X is a one-dimensional 
array containing the values of the vectors of variables U, V, and W (stored 
contiguously), T is the fire time, ¢, and Z is an array which receives and stores 
(also contiguously) the ouputs of F', G, H, and J evaluated at the input values of 
X and T. Upon receiving a call from the numerics, the communications section of 
the interface calls the physical subroutines to evaluate the components of Z. The 
communications section interprets between the numerics’ abstract representation 
of functions and variables as nameless components of vectors, and the physics’ 
representation of functions and variables as named entities, carrying physical 
significance. This ensures that the appropriate physical subroutine is called with 
the appropriate input variables to evaluate each component of Z. 


4.2.1 Subroutine UNPACK 


UNPACK translates a vector of variables stored in the format used by the 
numerics into the format used by the physics. It copies the values from a one- 
dimensional input array X into common /VAR/ so that the value in X(I) is 
assigned to the position of /VAR/ corresponding to the I*® variable in the system. 
Thus, if the radius of the fire on object 2 in room 3 is the 25'" variable in the 
system, i.e. the 25°» component of the vector (U,V, W, Y), then a call to UNPACK 
copies the value of X(25) into the element ZRFZZ(2,3) in /VAR/. 


The information stored in the integer array INSYS in common /SYSPTR/ gives 
UNPACK the correspondence between the variables in the system and positions 
in /VAR/. The value in INSYS(I) is an integer giving the location in /VAR/ 
‘corresponding to the I*® variable in the system. Subroutine SETSYS, in the 
control section of the interface, sets up the pointers in this array whenever the 
system of equations changes. 


4.2.2 Subroutine FUNCT 


Subroutine FUNCT is the numerics’ window to the interface. It takes T 
containing time and a one-dimensional array X containing values for the variables 
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in the system as input, and fills the one-dimensional output array Z with the 
values of the functions, associated with the equations in the system. FUNCT uses 
EVALFP and ERFDE to call the physical subroutines. 


FUNCT performs its task in three steps. First, a call to UNPACK loads the 
variable values from X into /VAR/ so that they may be passed to the physics as 
named, semantically meaningful entities rather than as nameless components 
of a vector. Second, FUNCT calls EVALFP and ERFDE which in turn call the 
physical subroutines to evaluate the required functions leaving the results in 
/VAR/. Finally, FUNCT copies the results of this function evaluation from /VAR/ 
into the output array Z, in the appropriate order. 


After a call to UNPACK, every non-eliminated variable in /VAR/ relevant to 
the current problem contains a meaningful value. The positions in /VAR/ cor- 
responding to non-eliminated variables which are components of X (variables sent 
to the numerics) receive values in the call to UNPACK. The positions in /VAR/ 
which correspond to variables that are relevant to the current problem but not 
varying, and therefore not sent to the numerics, retain fixed values assigned 
in SETSYS. For example, if an object is not burning, SETSYS sets the element 
of /VAR/ which corresponds to the pyrolysis rate TMOZZ for that object to 0.0 
so that subroutines which take TMOZZ as input return meaningful values. Note 
that positions in /VAR/ which correspond to variables not relevant to the cur- 
rent problem may contain meaningless values. For example, variables describing 
processes for objects which do not exist in the current run may have no assigned 
value. 


Subroutines EVALFP and ERFDE evaluate the necessary functions using inputs 
from /VAR/, and store the results in /VAR/. EVALFP calls all the physical 
subroutines associated with fixed-point equations while ERFDE calls those physical 
subroutines associated with root-finding and differential equations. As some of 
the subroutines called by ERFDE may take eliminated variables as input, the call 
to EVALFP must come before the call to ERFDE. (See the discussion of variable 
elimination in the section on ERFDE and EVALFP for a fuller explanation.) 


After the calls to EVALFP and ERFDE, /VAR/ contains the results of all the 
function calls. FUNCT packs these results into the output array Z using the 
information from the integer array PACK in common /PCKPTR/. The integer 
stored in PACK(I) is a pointer to the location in /VAR/ which holds the output 
of the function corresponding to the I*® component of Z. Thus, FUNCT plucks 
function outputs from /VAR/, accessed via the array XIN, and places them 
successively in Z. 
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4.2.3 Subroutines EVALFP and ERFDE 


Subroutines EVALFP and ERFDE call the physical subroutines to evaluate the 
functions associated with the system of equations. The first of these subroutines 
EVALuates the functions associated with the Fixed Point equations and the second 
Evaluates the functions associated with the Root-Finding and Differential Equations. 
The subroutines pass the parameters in common /BLDNG/ to the physical sub- 
routines and use the values in common /VARYNG/ to decide which physical sub- 
routines to call. Both subroutines take input values for the variables from com- 
mon /VAR/ and store the output function values back in common /VAR/. The 
flow of information in and out of /VAR/ is organized in such a way that certain 
physical subroutines receive as input the outputs of other physical subroutines. 
By selecting the appropriate calling sequence for the physical subroutines, the 
programmer may take advantage of this “composition” effect to reduce the size 
and/or complexity of the problem which the numerics must solve. This is 
described in the following section. 


4.2.3.1 Variable Elimination 


The cost of finding a numerical solution to a system of equations depends 
very heavily on the number of equations which must be solved simultaneously. . 
Analytically, one may be able to eliminate certain variables from a system 
of equations, thereby reducing the expense involved in solving that system. 
Using an analogous technique, the variable elimination scheme used in CFC 
VI takes advantage of the sparseness of the Jacobian matrix corresponding 
to the model’s equations to greatly reduce the size of the problem which the 
numerics must solve. This section gives a brief descritpion of the idea behind 
_ variable elimination and discusses in greater detail the current program’s variable 
elimination scheme and possible extensions for future programs. Note that a 
complete understanding of the function of the communications section of the 
interface requires a more thorough understanding of the concepts involved in 
variable elimination than this section attempts to provide. For a more detailed 
discussion of the theory of variable elimination see [1], [6], and [5]. 


Variable elimination applies only to fixed-point algebraic equations which 
are to be solved numerically. An example taken from [6] illustrates most of the 
important concepts. Suppose one wishes to solve the system 
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t1 = 9i(21, 22) 
Lo = go(23, Zs) 
3 = 93(21, 24) (4.8) 
L4 = ga(Zo, 25) 
ts = gs(21, 22) 


via Newton’s method. Most of the expense of this procedure occurs in solving 

the linear system J Az = ¢€ for Az, where J is the Jacobian matrix of the system 

of equations | 
F(21,%2,%3,%4,%5) = (%1 — 91(%1, 22), 22 — go(%3,25),..., 25 — 9s(21, 22), 


Az is the correction to be applied to the variables, and ¢ is the observed vector 
of divergences from the solution. In this case, J has the form 


x *x-0° 00 
Ol XU ex 
re] i Leese 
pa (MRO) = xO1x0 
J j= 1,...,5 OxOdOlx 
x x 0-0] 


where the x’s represent the (possibly) non-zero entries. Now instead of (4.8), one 
could solve the system 


21 = 93(21, 22) 
where g}(1, 72) = 91(21, 22) 
Za = 95(21, 22) 
where go(71, to) = g2(gs(21, 94(22, 95(21,22))), go(1,22)) 


23 > 93(21, 22) 


where 93(21, 22) = gs( 21, 94(2, 95(%1, t2))) ae) 


£4 >= 94(21, 22) 

where 93(21,22) = ga(z2, 95(21, Z2)) 
5 = 95(Z1, 22) 

where g5(21,22) = 95(21, 22). 


Here, the 9; functions are obtained from the g,; functions by substituting the 
right hand side of equations from (4.8) for the corresponding variables until only 
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Z; and Z2 appear on the right hand side of the resulting system. The solutions 
to (4.9) are exactly the same as those for (4.8). Moreover, the Jacobian matrix 
of the new system has the simpler form 


x x 
rex 

M Aga) Ba << ; (4.10) 
x x 
xx001l 

which allows the numerics to obtain values for Az at far less expense. The 

numerics need only solve the simultaneous equations 7; = 9:(21, Zo) and ty = 

9o(21, 22), after which values for z3, 24, and z5 may be deduced by substituting 

the known values for z; and zo in the remaining equations. 


In general, after variable elimination is applied to a system of equations the 
variables which do not appear on the right hand side of the resulting system are 
called eliminated or dependent variables. Those which remain on the right hand 
side are called non-eliminated or independent variables. The goal of any variable 
elimination scheme is to minimize the number of independent variables. 


Though the construction of the g; may seem complex as written above, in 
practice it is quite easy to construct a procedure which evaluates the 9; given 
the g;. In the above example, suppose Gi, G2,..., G5 are physical subroutines 
which evaluate gi,...,g5 respectively. Subroutine Gi takes one input for each 
input of g; and returns a single output corresponding to the value of g; at 
the input point. If V1, V2,..., V5 are FORTRAN variables then the follow- 
ing sequence of FORTRAN statements evaluates Gr yh ge given x; and Zo: 


C Start with x# stored in Vl and x# stored in V2. 

C The first two arguments of each Gi are the inputs, 
C the third is the output. 

c 


CALL G5(V1,V2,V5) 
CALL G4(V2,V5,V4) 
CALL G3(V1,V4,V3) 
CALL G1(V1,V2,V1) 
CALL G2(V3,V5,V2) 


C ay Fd 
C Vi1,....,V5 now contain 3, (Xs sXgoeeees de (x, 5%) 
C respectively. 


Figure 9. Variable elimination in FORTRAN. 
Thus, if the physical subroutines return their outputs to the same block of 
memory from which they take their inputs, one need only choose the correct call- 
ing order for the physical subroutines in order to implement a variable elimina- 
tion scheme. 


In the current program, the physical subroutines place their outputs in the 
same section of common /VAR/ from which they take their inputs, just as the 
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hypothetical physical subroutines in the above example use the same storage 
locations for input and output. Thus, by choosing the correct order of evaluation 
for the physical subroutines, the programmer can reduce the number of equations 
which the numerics must solve simultaneously. This calling order is stored in 
the array EVORD of subroutine EVALFP. 


Finding an order of evaluation which gives a small number of independent 
variables poses a formidable challange in large problems. In general, it is not 
possible to find the optimal evaluation order without examining every possible 
order. J. D. Ramsdell has written a program, JUGGLE, which suggests an evalua- 
tion order yielding a relatively small number of independent variables for most 
problems [5]. 


There are a number of ways in which variable elimination might be applied to 
the fire model, ranging from complex implementations which should produce the 
most time efficient results to simple ones which sacrifice some computation time 
but, nonetheless, are a vast improvement over the original system of equations. 


Only a certain subset of the set of all the model’s possible equations applies to 
the situation occurring at each time during a particular calculation. Moreover, 
at any given time in a calculation, the value of a particular physical function may 
depend on only a few of the variables which it takes as inputs. For example, 
in the current model, the flux from a source to a target object depends on 
a host of variables when neither object is immersed in the hot layer but is 
fixed at 0 (and therefore depends on no variables) when either of the objects 
is shielded by the opaque soot of the hot layer. The most efficient and most 
complex implementation of variable elimination would use all of this information 
to select the optimal evaluation order for the functions at frequent intervals. 
However, Ramsdell’s algorithm for finding an efficient evaluation order for the 
functions requires extensive manipulation of directed graph data structures, and 
uses recursive function calls and dynamic memory allocation. None of these 
features is readily mimicked in a FORTRAN program. 


Therefore, the current program does not actually comprise a dynamic im- 
plementation of Ramsdell’s algorithm for variable elimination. It cannot take 
advantage of the minute to minute changes in the interdependencies among the 
variables. Instead, it uses an evaluation order, chosen a priori using the separate 
program JUGGLE, based on the most general set of interdependencies which could 
ever occur among the variables in the model. As an example, consider the cal- 
culation of the flow rates through vents between adjoining rooms calculated in 
subroutine FLOW. If the interface between the hot and cold layers in both of the 
rooms is above the transom of the vent, then the flow rates through the vent 
depend only on the difference in the pressures measured at the floors of each of 
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the rooms. If, on the other hand, the layer interface in either room is below the 
transom of the vent, then the temperature and depth of the hot layer in that 
room affect the flow rates as well. But the current program does not take ad- 
vantage of the smaller set of interdependencies prevailing when the hot layer has 
not yet descended to the level of the vent. Instead, JUGGLE selects an evaluation 
order for the functions based on the assumption that the flows depend on the 
layer properties at all times.” 


Intermediate schemes which reset the order of evaluation more frequently 
have been suggested. One could run JUGGLE to get a function evaluation order 
for each of the several special cases that occur during a run. For example, one 
order can be used when a room contains no burning objects, a second order when 
it contains one burning object, a third when it contains two burning objects, 
and so on. At the moment, however, it appears that such schemes would result 
in very little time savings. In the future, as more sophisticated programming 
languages become more widely available, it should be possible to incorporate the 
full power of Ramsdell’s variable elimination scheme into the interface so that 
the function evaluation order is updated more frequently. 


Applying Ramsdell’s variable elimination procedure to the fire model equa- 
tions raises several difficulties. First, Ramsdell’s procedure properly applies only 
to systems of fixed-point equations whereas the fire model’s systems include root- 
finding and differential equations. To solve this problem, the program treats each 
of the variables corresponding to one of these root-finding or differential equa- 
tions as independent variables so that variable elimination applies only to the 
variables corresponding to fixed-point equations. Since, in fact, most of the equa- 
tions are fixed-point, algebraic equations, little is lost by this artifice. Moreover, 
the equation corresponding to any independent variable may be evaluated after 
all of the equations corresponding to dependent variables without upsetting 
the variable elimination scheme (regardless of how JUGGLE orders these evalua- 
tions with those of the dependent variables). Hence all of the root-finding and 
differential equations can be (and are) evaluated as a block by subroutine ERFDE, 
with the call to ERFDE occurring after the call to EVALFP. 


Second, some physical subroutines return more than one output per call. 
Such a subroutine must be called at the earliest point in the calculation at 
which any of its outputs is needed. Since such subroutines restrict the possible 
evaluation orders for the variables, in general they tend to increase the number 


2Using Ramsdell’s graph-theoretic terminology, the evaluation order for the functions is chosen 
to minimize the resulting size of the feedback vertex set for a graph which has the property 
that it contains as a subgraph the Newton graph of every system of equations ever occurring 
in an actual run. 
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COMMON /VAR/ TELZR(MR),TELZD(MR),ZMLZZ( MR), TMLZZ(MR),ZELZZ(MR), 
TELZZ(MR),ZHLZZ(MR),ZKLZZ(MR),ZYLOZ(MR),ZYLDZ(MR),ZYLMZ(MR), 
ZYLSZ(MR),ZYLWZ(MR),ZPRZZ(MR),ZULZZ( MR), TEPSM(MR),TEOSM(MR), 
TMODS( MR), TMOMS(MR), TMOSS( MR), TMOWS(MR),ZMDL(MR), 
ZMML(MR),ZMOL(MR),ZMWLOMR),ZMSL(MR),AVLSM(MR),WVLS(MR), 
PPLR(MR),TMOSM(MR),TMPSM(MR),TMUSM(MR),TMDSM(MR),TEUSM(MR), 
EPLSM(MR),BOXYS(MR), TMDZF (MR), TMWZF(MR),TMSZF(MR), TMMZFCMR), 
TMOZF(MR),FQPOR(MO,MR),ZKOZZ(MO,MR),ZMOZZ(MO,MR), 
TMOZZ(MO,MR),TEOZZ(MO,MR),ZHPZZ(MO,MR),TMPZZ(MO,MR), 
TEPZR(MO,MR),ZRFZZ(MO,MR),ZHFZZ(MO,MR),FQLOR(MO,MR), 
FQWOR(MO,MR),TMOZB(MO,MR),PHI(MO,MR),AQWZZ( MW) ,ZKWZZ(2,MW), 
FQLWR(2,MW),FQPWR(2,MW),FQLWD(2,MW), 
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1 TMNET(MR), 
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1 DZMLZZ(MR),DZELZZ2(MR),DZ2MOZZ(MO,MR),DZ2MDL(MR),DZMML(MR), 
2 DZMOL(MR),DZMWL(MR),DZMSL(MR),DAQWZZ( MW) ,DZRFZZ(MO,MR) 


Figure 10. Common /VAR/. 


+eWON FATA WN = 


of independent variables. Such subroutines should therefore be avoided where 
possible. 


Since the current program uses a fixed calling order for the physical sub- 
routines, each fixed-point equation is either an eliminated variable equation or a 
non-eliminated variable equation independent of which equations and variables 
happen to be relevant to a particular calculation. Thus, each equation is of one 
of four types: (1) root-finding equation; (2) differential equation; (3) fixed-point 
equation for non-eliminated variable; (4) fixed-point equation for eliminated vari- 
able. Therefore, when SETSYS determines that a particular equation is in the 
system, there is never any ambiguity as to what type of equation it is. (In the 
future, if the subroutine calling order is updated periodically in the course of 
a run, SETSYS will have to place a call to some subroutine performing the func- 
tion of JUGGLE at each discontinuity in order to determine which of the relevant 
fixed-point equations are eliminated variable equations in the context of the new 
system.) 


4.2.3.2 Organization of the Subroutines 


Subroutines EVALFP and ERFDE take input variables from and return output 
function values to common /VAR/. The statement describing /VAR/ for EVALFP 
and ERFDE appears in figure 10. 


All of the arrays by which EVALFP and ERFDE access /VAR/ are of the REAL*8 
data-type. This block of REAL*8 memory locations is conceptually divided into 
three sections, delimited in figure 10 by the comment lines of “x”s. When 
FUNCT first calls EVALFP, section 1 of /VAR/ contains a meaningful value for 
each of the variables relevant to the current problem. Sections 2 and 3 of 
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/VAR/ contain garbage. Subroutine EVALFP calls each of the physical subroutines 
corresponding to fixed-point equations relevant to the current problem. It passes 
variable values from section 1 of /VAR/ and places the results back in section 1 
of /VAR/ immediately upon exiting from the subroutine. These new values may 
then be passed along to other physical subroutines, thereby facilitating variable 
elimination. 


After EVALFP completes its evaluation of the fixed-point equations, control 
passes to ERFDE which evaluates the functions associated with the root-finding 
and differential equations. Like EVALFP, ERFDE passes input variables from 
section 1 of common /VAR/ to the physical subroutines. However, the outputs 
of the subroutines corresponding to differential equations go into section 3 of 
/VAR/ while the outputs of subroutines corresponding to root-finding equations 
go into section 2 of /VAR/. (Section 2 of /VAR/ contains one memory location 
earmarked to receive the output of each root-finding equation. In the current 
program, TMNET(R) receives the output of the function associated with the mass 
flow equation for room R. Section 3 of /VAR/ contains one memory location 
earmarked to receive the output of the function associated with each differential 
equation. For example, DZMML(R) receives the time rate of increase of the mass 
_ of carbon monoxide in the hot layer of room R.) Since no physical subroutine 
receives the contents of sections 2 and 3 of /VAR/ as input, the order in which 
ERFDE calls the physical subroutines is unimportant. 


At a given time in the run, EVALFP and ERFDE need not call every possible 
physical subroutine. Calling subroutines which describe processes not relevant 
to the current problem is a waste of time. Still worse, such a subroutine call may 
cause an error. For example, a call to a subroutine which attempts to calculate 
the flame velocity on a cold object will, at best, return garbage, and at worst may 
present the subroutine with a set of values which it cannot handle. Therefore, 
each subroutine call in EVALFP and ERFDE is predicated on the information which 
SETSYS places in common /VARYNG/. The 78 position in common /VARYNG/ 
contains a 0 if the i*" variable in /VAR/ has a fixed value at the current stage 
in the calculation or has no meaningful value, and contains an integer n > 0 
if the n*® possible physical equation corresponds to that variable in the current 
problem. Thus, EVALFP and ERFDE call a physical (hypothetical) subroutine, 
SUBR, if and only if SUBR computes the function associated with one of the 
currently relevant equations. 
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Appendix A: 
Introducing a New Physical Subroutine 


CFC VI is designed to facilitate frequent changes in the set of physical 
equations and variables. This section serves as a practical guide for introducing 
new equations and variables into the model. The discussion assumes a general 
familiarity with the concepts and terminology introduced in the main body of 
the document. 


Modifying the Existing Equations 


The programmer may make any desired changes to the existing set of physical 
subroutines with the following considerations: 


1. If the resulting subroutine depends on the same physical parameters and 
variables as the old subroutine and applies to the same physical situations, 
no further modifications are necessary. 


2. Ifthe resulting physical subroutine depends on additional physical parameters, 
the programmer must update the argument list of the subroutine and the cor- 
responding argument list in the call to the subroutine in EVALFP and ERFDE. 
No other modifications need be made. 


3. If the resulting physical subroutine has a different range of applicablity 
than the old subroutine, SETSYS and NWSTAT must be updated to reflect 
the change. For example, if the programmer adds a new case to TMPO2 so 
that it can compute the surface temperature of a burning object as well as 
that of a heating object, SETSYS must be adjusted so that it no longer fixes 
the surface temperature of an object at the ignition temperature when the 
object ignites. ree 


4. If the resulting physical subroutine depends on more or serase he physical 
variables than the old, the variable elimination scheme should be updated. 
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This is properly accomplished with the aid of JUGGLE. However, if JUGGLE 
is not available or if the programmer is in a hurry, (s)he may temporarily 
change the values in array EVORD of subroutine EVALFP so that the affected 
subroutine is called last, and update SETSYS so that it does not treat any 
of the outputs of the affected subroutine as eliminated variables. If the 
programmer is relatively fluent in the concepts involved in variable elimina- 
tion, other procedures may be apparent. For example, if the variables added 
to the calling list are not eliminated variables, no changes need be made in 
the subroutine calling order. 


Adding a New Physical Subroutine 


If the programmer wishes to introduce a new physical subroutine into the 
model, (s)he must update EVALFP or ERFDE to call the subroutine under the 
appropriate circumstances, and modify SETSYS to specify those circumstances. 
For example, if the programmer wishes to introduce a new physical subroutine 
to compute the surface temperature of a burning object, (s)he must add a call 
for this subroutine in the section of EVALFP which evaluates ZKOZZ and predicate 
this call on the value of FZKOZZ in common /VARYNG/; for example EVALFP might 
call the new subroutine for object 0 in room R only if FZKOZZ(0,R) .EQ.2 and 
must call the old subroutine only if FZKOZZ(0,R) .EQ.1. In addition, subroutine 
SETSYS must be updated so that it places a 1 in the appropriate postion of 
common /VARYNG/ when the object is heating, and a 2 in that position when the 
object is burning. 


Whenever a new subroutine computing output variable V takes input vari- 
ables which previously did not affect the value ofV, the variable elimination 
scheme must be updated as described in “Modifying the Existing Equations,” 
above. 


Adding a New Physical Variable 


Whenever the programmer introduces a new variable into the model, (s)he 
must make the following changes: 


1. Add a REAL*8 array to section 1 of common /VAR/ containing one element 
for each entity to which the physical variable applies. The name for this ar- 
ray should not exceed five characters. For example, if one wishes to calculate 
the temperature of the floor in each room, a REAL*8 array ZKFZZ(MAXRMS) 
must be added to /VAR/, where ZKFZZ(R) corresponds to the floor tempera- 
ture in room R and MAXRMS is the maximum number of rooms allowed in 


4, 


T. 
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the model. Section 1 of /VAR/ is organized so that all of the arrays indexed 
XXXXX(MAXRMS) appear first, all those indexed XXXXX(MAXOBS , MAXRM8) ap- 
pear next, all those indexed XXXXX(MAXWLS) appear next, and all those in- 
dexed XXXXX(2,MAXWLS) appear last. The programmer should place the new 
array with the others having the same type of index. 


Increment the parameter (RMVARS, OBVARS, WLVARS, or SDVARS) which gives 
the number of arrays in /VAR/ having the same indexing as the new variable. 
These parameters are currently assigned in the include file DIMN.INC. 


Insert a five character string giving the name of the new array in the ap- 
propriate place in the DATA statement initializing NAMES in the CHARACTER 
FUNCTION VARNAM. The names in this statement are in the same order as the 
arrays in common /VAR/. 


In order to output values for the variable, modify WRIT to write out the 
variable (after computing it via a call to the physics if it is not one of the 
variables included in the array X). 


Add a new REAL*8 array to common /SCALER/ in BLOCK DATA SCALER. 
This array should have the same name as that introduced into /VAR/ but 
suffixed with an “A”. The arrays in /SCALER/ are in the same order as their 
counterparts in /VAR/. Add a new DATA statement to BLOCK DATA SCALER to 
initialize the new array with a value giving the smallest physically significant 
value for that variable. 


Add an integer array to common /VARYNG/ in subroutines ERFDE and EVALFP 
corresponding in size and position to the new array in /VAR/. This array 
should have the same name as that introduced into /VAR/ but prefixed with 
an “eo : 


Add the new variable to the argument lists in EVALFP and ERFDE for any 
subroutines using it as input or output. 


Udpate NWSTAT so that the new variable gets a proper initial value at each 
discontinuity. 


Introduce a new integer parameter into the parameter scheme in SETSYS, 
which gives the postion in /VAR/ of the first word of the new array. The new 
parameter should have the same name as the array in /VAR/ but prefixed 
with a “P”. For example, having introduced the array ZKFZZ(MAXRMS) into 
/VAR/ between TMOZF and FQPOR, one should change the block of PARAMETER 
statements in SETSYS from 
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12. 


13. 


14. 


PTMMZF 
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= PIMSZF + MAXRMS 
PTMOZF = PTMMZF + MAXRMS 
PFQPOR = PIMOZF + MAXRMS 
PZKOZZ = PFQPOR + MAXRMS*MAXOBS 

to 

PTIMMZF = PTMSZF + MAXRMS 
PTMOZF = PIMMZF + MAXRMS 
PZKFZZ = PTMOZF + MAXRMS 
PFQPOR = PZKFZZ + MAXRMS 

= PFQPOR + MAXRMS*MAXOBS 


PZKOZZ 


Update SETSYS so that the new variable is active and sent to the numerics 
at the proper times. 


If a differential equation ever describes the progress of the new variable, add 
a corresponding array to section 3 of /VAR/ reserved to receive the outputs 
of the function associated with this equation. To get the name for this 
array, prefix the name of the new array in /VAR/ with “D”. Similarly, if a 
root-finding equation ever describes the progress of the new variable, add a 
corresponding array to section 2 of /VAR/ reserved to receive the output of 
the function associated with this equation. The name for this array should 
be descriptive of the output of the function, e.g. TMNET for Time derivative 
of NET Mass in a room. 


In case (11) applies, update the parameter block in SETSYS, in a manner 
analogous to that described in (9), to include a pointer to the new array. 


Update the subroutine calling order, fixed by the DATA statement for array 
EVORD in subroutine EVALFP. This is properly accomplished using JUGGLE. If 
JUGGLE is not available, evaluate the new variable last and do not attempt 
to make it an eliminated variable. In this case, the program will run, though 
probably not in the most efficient possible way. 


Recompile all modules. 
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Appendix B: 
Numerical Considerations for Fire Models 


by John D. Ramsdell 


Introduction 


In any mathematical model which one expects to solve numerically, one 
must make certain restrictions on the form of the equations in the interest of 
developing a model that can yield useful answers using reasonable amounts of 
computer resources. This note will give some guidelines, based on continuation 
methods, for producing useful mathematical models. 


A fire model will be thought of as a system of m differential equations and n 
algebraic equations 


(1) = Hal), (,8), (ta) = 2°; 


0 = (x(t), y(?), 2); (1) 


where z:(R-R™ and y: R-R”. Usually f and g do not depend explicitly on 
ic 


It is assumed that equation (1) has a unique and continuous solution, and is 
well-posed with respect to initial conditions. Well-posed means small perturba- 
tions in the initial conditions result in small changes in the solution. 
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Equation (1) can be viewed as the following continuation problem: 


bar ls t) — (te fhe f(z, y, t)dt — 2{to) 


9(z, y, t) 
(2) 


= 0. 


Equation (2) is put in a more standard form by defining z to be the con- 
catenation of z and y, so equation (2) becomes F(z,t) = 0. The Jacobian, 
J(z,t), of F(z, t) is defined by 


OF; 
Jig(2,t) = *(z, 1). 
x) 


The concept of a norm will be used extensively in the remainder of this note. 
Vector norms have properties intuitively associated with the idea of length. Any 
vector norm satisfies 


|z|| > 0 iff z 4 0; 
jaz] = |a|||z/], a a scalar; 
Iz + yllS[l2l] + [ly 


A common norm is the 2-norm which is the square root of the sum of the 
squares of each component of the vector: 


3 
lalla = (28) - 


A norm often used in practice is the infinity norm which is the absolute value of 
the component which has the largest magnitude: 


I2Ilo0 = max|= 


A matrix norm can be defined by a vector norm. For a matrix, A, 


||A]| = max Aa 
20 ||2| 
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This norm has the following properties: 


||Aj| > 0, iff A 5 0; 
||a.A}| == |a| |LAl|, a a scalar; 
A+ Bis lAll + Bil; 

|| ABI] <All BI; 

||Az||< IAI] [l2I. 


The 2-norm of a matrix is the square root of the largest eigenvalue of AYA, 
and the infinity norm is 


Ax 
|Alloo va wt a oa max ) Ag 


The Conditions 


Numerical methods for solving continuation problems often assume the fol- 
lowing three conditions on F(z, t) [7]. 


I. <A “good” initial guess is available at each time step. 
I. For some constants kK, L, and M; z(t) a solution of F(z,t) = 0; and 7 near 


a(t), 


B  |lJ(2,t)— J(Z,t)I| < Liz — al 
CG |ju(z,t)Ih |] I-*(2, | = cond( J(z, t)) < M. 


fit. Condition If holds in a sufficiently large ball around the solution. 


Condition I guarantees that a global search procedure is not needed to find a 
solution. A good guess can be found in the fire model by making the time step 
small enough, assuming equation (2) has a continuous solution. 


Condition I! bounds the size of some important quantities. Differentiating 
equation (2) with respect to time gives a system of differential equations useful 
in explaining the bounds: 


=() = == JiT*(z, 1) (2, t). (3) 
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Condition IIA bounds the size of (OF'/At)(z, t). This contributes to bounding 
the size of the derivative of z. Very small time steps may be required to solve 
equation (2) when this bound is not satisfied. This result and the following ones 
are independent of the particular numerical method being used. 


Condition IIB requires that the Jacobian is Lipschitz continuous, in other 
words, the Jacobian is not too sensitive to changes in z. This is important when 
finite differences are used to estimate the Jacobian. It justifies the use of large 
enough increments in z to avoid errors due to round-off. Conditions IIA and IIB 
combined at times justify using the Newton chord method to solve equation (2). 
This method uses the same estimate for the Jacobian at many time steps. 


Condition IIC bounds the condition number of the Jacobian. The condition 
number measures the sensitivity of the time derivative of z to errors in the 
Jacobian. This can best be explained by looking at the following linear system: 


Az = b. 
Assume A is perturbed so that instead we solve 
AU + F\(z+ Ar) = 6, 
Haw \|F'|| < 1. We would like to relate changes in A (via F’) to changes in z 
(via Az). 
Ar =(A(U+ F))~'b— A716 
=([1+F)'A—'b— Ab = (7 +F)y*— I)z 
=(+F)70—-(U4+ F))t=—-U4+ F)-'F 2 


thus 


vay SU +P) 


Since F = A~!AA (where AA = A(J + F) — A), 


|AAl 
[|All 


[AAI _ 


| 
ond(A 
Al 4 'TAL 


FI < ATI [Al] = (Ila jaa? HAs 


It can be shown that 


1 
P-L Pier lee 
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so we get the desired result of 


ai]. _ cond(A)'rayh re 
{Ele es cond(A) inal 


if 
|AAll 


conde baaip ile 


Equation (4) shows that the condition number of a matrix measures the 
sensitivity of the solution to changes in the matrix. Returning to equation (3), 
we can write it in a form like a linear system: 


da — oe t). 


J(2, 1) At 


Bounding the condition number of J makes sure that the derivative of z is not 
too sensitive to variations in J. 


Condition II allows a program to assume condition II even at some 2’s that 
do not solve equation (2). In particular, condition I is satisfied on the path from 
the initial condition to the solution at each time step. What we really need is 
that the conditons be satisfied just at the places we evaluate Ff. Thus, if we 
choose one-sided differences to estimate J, we do not need condition II to hold 
in the other direction. 


Application to CFC VI 


The current fire model is being solved using an algorithm for algebraic- 
differential systems developed by Gear [2,3]. When the fire model satisfies 
conditions I-III, Gear’s code quickly gives a solution to the model. However, 
at. some points of the model, the code runs poorly. We know the Jacobian is 
discontinuous at some points, but we are not sure what happens at other points. 
It is important to find out if the conditions are being violated at these points; 
for if they are, we cannot expect any numerical procedure to find a solution, and 
we must change the physics. If the conditions are satisfied, we should focus on 
finding new numerical algorithms to solve equation (2). 


It is hoped that future writers of fire models will strive to create models that 
satisfy the conditions as much as possible. Good solutions can be expected from 
these models. 
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An Alternative Solution Method? 


Some people suggest equation (1) is best solved by converting it into a system 
of ordinary differential equations. Differentiating g(z, y,¢) with respect to ¢ gives 


dz 
— = t 


Lh a (2 


, 99 
dt OY ot +. 22 7(2,4,1)), 


Experience has shown that these equations are usually harder to solve than 
directly solving equation (1). 
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Appendix C: 
Dictionary 


The following is an alphabetical list of the most important terms and symbolic 
names associated with CFC VI, along with a little information on each. Included 
in the list are SUBROUTINE, FUNCTION, PROGRAM, ENTRY, and BLOCK DATA names; 
COMMON block names; the names of arrays contained in COMMON block areas; 
PARAMETERS which are used throughout the program; the names of the most 
important local variables; the names of other programs used in conjunction with 
CFC VI; and terms used frequently in the documentation. 


ABSRB1i Calculates the infra-red absorptivity of the hot layer. 
SUBROUTINE in the physics module. 


ABSRB2 Same as ABSRB1 but more sophisticated. 
SUBROUTINE in the physics module. 


ABSRB3 Same as ABSRB2 but still more sophisticated. 
SUBROUTINE in the physics module. 


ADINIT (Algebraic and Differential equation solver INITialization.) Called 
whenever the system of equations changes to restart the in- 
tegrator. 

ENTRY to ALDIFF. 


AFIRE(0,R) Fire spread parameter for the object. 
REAL*8 in COMMON /BLDNG/. 
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ALDIFF (ALgebraic and DIFFerential equation solver.) Solves the sys- 
tem of algebraic and differential equations. See chapter 3. 
SUBROUTINE in numerics module. 


algebraic equa- ; 
tion Any equation of the system which is not a differential equation. 
See sections 1.2.1 and 2.1.2. 


ALPHA Plume entrainment coefficient. 
REAL*8 in COMMON /BLDNG/. 


ANGH (0,8) Angle of the surface of an object with the horizontal plane. 
REAL*8 in COMMON /BLDNG/. 


ANGV(0,R) Angle of the surface of a vertical object with the vertical plane 


through the x-axis; i.e., with the X-Z plane 
REAL*8 in COMMON /BLDNG/. 


AQWZZ(W) Heat which has been absorbed by a unit area of the wall 
material of wall W. 
REAL*8 in section 1 of COMMON /VAR/. 


ASYMP Used by ABSRB3. 
SUBROUTINE in the physics module. 


AVLSM (R) Total area of the vents in the walls of room R which is wetted by 
the hot layer. 
REAL*8 in section 1 of COMMON /VAR/. 


BETA Unused in CFC VI 
REAL*8 in COMMON /BLDNG/. 
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BLDNG 


BOXYS (R) 


BRNOUT 


BURN 


BURNER 


CFC V 


CHEBY 


CHEM 


CHI(O,R) 
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(Describes the BuiLDING to be burned.) Common block which 
contains all the values needed by the physics module but not 
treated as variables by the numerics module. For example, 
this common block contains the dimensions of rooms in the 
building and the positions of objects in those rooms. See 
introduction to chapter 4. 

COMMON block. 


(Burned OXYgen Sum.) Total rate of burning of oxygen in the 
plumes of all the objects in room R. 
REAL*8 in section 1 of COMMON /VAR/. 


Value used by NWSTAT and UNPINT to indicate that an object 
has burnt out. 
INTEGER PARAMETER. 


Calculates the mass rate of burning of the fuel from an object. 
SUBROUTINE in the physics module. 


Used as a value for FIRTYP(0,R) to indicate that the object is 
a burner fire. 
INTEGER PARAMETER. 


Vent flow coefficient for all the vents. 
REAL*8 in COMMON /BLDNG/. 


The predecessor of CFC VI. 


Used by ABSRB3. 
SUBROUTINE in the physics module. 


Calculates the rate at which chemical energy is released from 
a burning object. 
SUBROUTINE in the physics module. 


Fraction of the heat of combustion released in the burning of 
the object. 
REAL*8 in COMMON /BLDNG/. 
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CHRCFR 


CNSUM 


CNTCIG 


CNVL1 


CNVW1 


COLD 


communications 
section of inter- 
face 


control/system 
set-up section of 
interface 


Used as a value for STATE(O,R) to indicate that the object is 
a charcoal fire. 
INTEGER PARAMETER. 


Used as a value for STATE(O,R) to indicate that the object has 
been consumed. 
INTEGER PARAMETER. 


Value used by NWSTAT and UNPINT to indicate that an object 
has ignited by contact with the flames of another object. 
INTEGER PARAMETER. 


Finds the rate of loss of energy from the hot layer in a room 
due to convective heating of the walls and ceiling. 
SUBROUTINE in the physics module. 


Calculates the convective heat transfer to one side of a wall 
from the adjacent hot layer. 
SUBROUTINE in the physics module. 


Used as a value for STATE(O,R) to indicate that the object is 
cold. 
INTEGER PARAMETER. 


The section of the interface which provides communications 
between the physics and the numerics modules. See sections 
4.2 and 1.2.4.2. 


That section of the interface module which controls the in- 
tegrator and the system of equations. Sometimes abbreviated 
to control section of the interface. See section 4.1 and 1.2.4.2. 
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control section 
of interface 


CORCON 


CORECT 


CP 


CVWID 


DAQWZZ (W) 


DECOMP 


dependent vari- 
able 


DERS 
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See control/system set-up section of interface. 


Can be used to write out a list of all those variables which fail 
to meet the corrector convergence criterion. 
SUBROUTINE. 


Performs the corrector cycle for the numerical package. At- 
tempts to iteratively solve the corrector equation given an 
initial guess and an estimate of the Jacobian matrix. See 
section 3.1.3. 

SUBROUTINE in numerics module. 


Specific heat of air. 
REAL*8 in COMMON /BLDNG/. 


Calculates the length of the line along which the hot layer 
intersects a vent. 
SUBROUTINE in the physics module. 


Time derivative of the total heat absorbed by a unit area of 
the material of wall W. 
REAL*8 in section 3 of COMMON /VAR/. 


Decomposes a real matrix by Gaussian elimination and es- 
timates the condition number of the matrix. 
SUBROUTINE in numerics module. 


An eliminated variable. One of the variables which is not in 
the system of equations which must be solved simultaneously 
by the numerics. See section 4.2.3.1. 


The index in the array X of the last integrated variable in the 
system. 
Local INTEGER in ALDIFF. 
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differential equa- 
tion 


discontinuity 


DLECK 


DSCTYP 


DSKOUT 


DXDCNC 


DXDPRD 


DZELZZ (R) 


DZMDL (R) 


DZMLZZ(R) 


An equation of the system with the form of equation (case b) 
(2.5). See section 2.1.2. 


An interruption which involves a change in the system of equa- 
tions to be solved and/or a discontinuous change in the values 
of the variables. See section 4.1.1. 


Used by ABSRB3. 
REAL*4 FUNCTION in the physics module. 


Stores the code for the type of discontinuity which UNPINT has 
detected. 
Local INTEGER in UNPINT, MARK6, and NWSTAT. 


The time between outputs to the disk file FIRE.DAT. 
Local REAL*8 in MARK6. 


Calculates the mass concentration of carbon dioxide in a hot 
layer. 
SUBROUTINE in the physics module. 


Calculates the rate at which the fire of an object produces 
carbon dioxide. 
SUBROUTINE in the physics module. 


Time derivative of the total energy of the hot layer in room R. 
REAL*8 in section 3 of COMMON /VAR/. 


Time derivative of the mass of carbon dioxide in the hot layer 
of room R. 
REAL*8 in section 3 of COMMON /VAR/. 


Time derivative of the mass of the hot layer in room R. 
REAL*8 in section 3 of COMMON /VAR/. 
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DZMML (R) 


DZMOL (R) 


DZMOZZ (0 ,R) 


DZMSL (R) 


DZMWL (R) 


DZRFZZ(0,R) 


EB(O,R) 


EGAS 


eliminated vari- 
able 
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Time derivative of the mass of carbon monoxide in the hot 
layer of room R. 
REAL*8 in section 3 of COMMON /VAR/. 


Time derivative of the mass of oxygen in the hot layer of room 
R. 
REAL*8 in section 3 of COMMON /VAR/. 


Time derivative of the remaining fuel mass of object 0 in room 
R. 
REAL*8 in section 3 of COMMON /VAR/. 


Time derivative of the mass of soot in the hot layer of room 
R. 
REAL*8 in section 3 of COMMON /VAR/. 


Time derivative of the mass of water vapor in the hot layer of 
room R. 
REAL*8 in section 3 of COMMON /VAR/. 


Time derivative of the radius of the fire on object 0 in room 
R. 
REAL*8 in section 3 of COMMON /VAR/. 


Emmisivity/absorptivity of the surface of the object. 
REAL*8 in COMMON /BLDNG/. 


Used by ABSRB3. 
REAL*4 FUNCTION in the physics module. 


A dependent variable. One of the variables which has been 
eliminated from the right hand side of the set of equations so 
that it does not appear in the system of equations which must 
be solved simultaneously by the numerics. See section 1.2.2 
and 4.2.3.1. 
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EMSVTY 


ENTPL 


EPLSM (R) 


EPS 


ERFDE 


ERROR 


EVALFP 


EVJAC 


_EXTING 


FCO(0,R) 


Used by ABSRB3. 
SUBROUTINE in the physics module. 


Calculates the rate of entrainment of layer gases into the upper 
part of the plume over an object. 
SUBROUTINE in the physics module. 


Total rate of entrainment of layer gases into the upper parts 
of all plumes in room R. 
REAL*8 in section 1 of COMMON /VAR/. 


The one step error tolerance for the integration. 
REAL*8 PARAMETER in MARK6. 


Evaluates the functions associated with the Root-Finding and 
Differential Equations in the system. See section 4.2.3. 
SUBROUTINE in communications section of interface. 


Used to store an estimate of the one-step truncation error for 
each of the variables in the system. 
Local REAL*8 array in ALDIFF and MARK6. 


EVALuates the functions associated with the Fixed Point equa- 
tions in the system. See section 4.2.3. 
SUBROUTINE in communications section of interface. 


(EValuate JACobian.) Uses finite differences to produce an 
estimate of the Jacobian matrix of the system of equations. 
See section 3.1.3. 

SUBROUTINE in numerics module. 


Used as a value for STATE(0,R) to indicate that the object has 
been extinguished. 
INTEGER PARAMETER. 


Mass fraction of carbon monoxide evolved in the pyrolysis 
gases of the object. 
REAL*8 in COMMON /BLDNG/. 
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FCO2(0,R) 


FH20(0,R) 


FINISH 


FIROOM(R) 


FIRTYP (0 ,R) 


fixed 


fixed-point equa- 
tion 


fixed-point vari- 
able 


FLAME 
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Mass fraction of carbon dioxide evolved in the pyrolysis gases 
of the object. 
REAL*8 in COMMON /BLDNG/. 


Mass fraction of water vapor evolved in the pyrolysis gases of 
the object. 
REAL*8 in COMMON /BLDNG/. 


The stopping time for the run. 
Local REAL*8 in MARK6. 


Indicates if there is a FIre in the ROOM R. 
LOGICAL in COMMON /BLDNG/. 


The type of fire on the object (GROWNG, POOL, or BURNER) if the 
object is burning. 
INTEGER in COMMON /BLDNG/. 


A physical variable is fixed if its value is not changing during 
a particular portion of the progress of a fire. See section 4.1.4. 


When discussing the physical subroutines, this term refers to 
any equation with the form of (case c) of equation (2.5). When 
discussing the numerics and interface, it refers only to those 
equations represented by the components of equation (1.3). 
(Notice that the components of (1.4) actually fit the form of 
(case c) of equation (2.5) as well. 


A variable whose progress is described by a fixed-point equa- 
tion. One of the components of W in equation (1.3). 


Used as a value for STATE(O,R) to indicate that the object is 
flaming. 
INTEGER PARAMETER. 
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FLOW 


FLUX 


FNUM 


FQLOR(0,R) 


FQLWD (8, W) 


-FQLWR(S,W) 


FQPOR(0,R) 


FQPWR(S,W) 


FQWOR (0,R) 


FS(0,R) 


Calculates the energy and mass flows through a vent. 
SUBROUTINE in the physics module. 


Calculates the flux incident on the surface of a burning target 
object. 
SUBROUTINE in the physics module. 


The number of root-finding equations in the system. 
Local INTEGER in the interface and numerics. 


Radiative flux from the hot layer in room R to object 0 in that 
room. 
REAL#8 in section 1 of COMMON /VAR/. 


Convective flux from the hot layer adjacent to side 8 of wall W 
to the surface of the wall. 
REAL*8 in section 1 of COMMON /VAR/. 


Radiative flux from the hot layer adjacent to side § of wall W 
to the surface of that wall. 
REAL*8 in section 1 of COMMON /VAR/. 


Flux from the plumes of all objects in room R to the surface 
of object 0 in room R. 
REAL*8 in section 1 of COMMON /VAR/. 


Radiative flux from all plumes visible to side 8 of wall W to the 
surface of the wall. 
REAL*8 in section 1 of COMMON /VAR/. 


Radiative flux from the walls in room R to object O in that 
room. 
REAL*8 in section 1 of COMMON /VAR/. 


Mass fraction of smoke evolved in the pyrolysis gases of the 
object. 
REAL*8 in COMMON /BLDNG/. 


v& 


FUNCT 
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Calls EVALFP and ERFDE in order to evaluate the FUNCTions 
associated with the equations in the system. This subroutine is 
called by the numerics and translates back and forth between 
the numeric’s data structures and those of the physics. See 
section 4.2.2. 

SUBROUTINE in communications section of interface. 


function associated 


with a physical 
equation 


functional sub- 
routine 


GNUM 


GROWNG 


HEATIG 


HEATNG 


HEIGHT 


The function ¢ in a physical equation represented in the form 
of equation (2.5). See section 2.2.1. 


A physical subroutine whose outputs depends only on its im- 
mediate inputs and not on any previous inputs. See section 
2:2 2k. 


The number of differential equations in the system. 
Local INTEGER in the interface and numerics. 


Used as a value for FIRTYP(0,R) to indicate that the object is 
burning as a growing fire. 
INTEGER PARAMETER. 


The size of a time step. 
Local REAL*8 in ALDIFF and MARK6. 


Value used by NWSTAT and UNPINT to indicate that an object 
has ignited due to heating of its surface. 
INTEGER PARAMETER. 


Integer parameter used as a value for STATE(O,R) to indicate 
that the object is heating. 
INTEGER PARAMETER. 


Calculates the height of the flames over an object. 
SUBROUTINE in the physics module. 
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HMIN The minimum allowable time step. 
REAL*8 PARAMETER in MARK6. 


HNUM The number of fixed-point equations in the system. 
Local INTEGER in the numerics and interface. 


HSTART The size of the time step which is used to restart the integra- 
tion after each call to ADINIT. 
PARAMETER in MARK6. 


HSTORE The step size at which the derivatives in X are scaled. 
Local REAL¥8 in ALDIFF and MARK6. 
independent vari- 
able A non-eliminated variable. One of the variables which is in the 


system of equations that must be solved simultaneously by the 
numerics. See section 4.2.3.1. 


INIT Reads data from the input file at the start of the run. 
SUBROUTINE. 
initial values This can refer either to the initial values needed to deter- 


mine a unique solution to a differential equation, or to the 
approximate first values for the variables which the numerics 
require in order to solve the system of equations. Even when 
initial values are provided for all the integrated variables, the 
numerics will not be able to solve the system unless a good 
first approximation to the solution is also provided. See sec- 
tion 4.1.3. 


input paramet- 
er In the context of a particular physical subroutine, a parameter 
which is input to that subroutine. 


input variable In the context of a particular physical subroutine, a variable 
which is input to that subroutine. 


&6 


integrated vari- 
able 


interface 


interface mod- 
ule 


interruption 


INUM 


IRTTY 


IWDSK 


IWTTY 


JAC 


JUGGLE 
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A variable whose progress is described by a differential equa- 
tion. One of the components of V in equation (1.2). 


See interface module. 


The section of the program which provides communication 
between the physics and the numerics modules, controls the 
integrator, and selects the appropriate system of equations to 
model a particular fire. 


An event which requires that the integrator converge at the 
time of its occurrence. See section 4.1.1. 


The number of eliminated variable equations in the system. 
Local INTEGER in the interface and numerics. 


Channel number for reads from the terminal screen. 
INTEGER PARAMETER. 


Channel number for writes to the disk file FIRE.DAT. 
INTEGER PARAMETER. 


Channel number for writes to the terminal screen. 
INTEGER PARAMETER. 


Storage area used to hold the Jacobian matrix of the system 
of equations. 
Local REAL*8 in ALDIFF, MARK6, and CORECT. 


PASCAL program written by J. D. Ramsdell which is useful in 
selecting the evaluation order for the variables so as to yield a 
relatively small number of independent variables. See section 
4.2.3.1 and [6], [5]. 
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kink 


LAYDPT 


LAYR (R) 


limit check 


LOOKUP 


MARK6 


MAXOBS 


MAXODR 


MAXRMS 


MAXVRS 


A discontinuity in a first partial derivative of one of the physi- 
cal equations. See section 2.1.3. 


Calculates the depth of the hot layer in a room. 
SUBROUTINE in the physics module. 


Indicates if there is a hot layer in the room. 
LOGICAL in COMMON /BLDNG/. 


A section found at the end of most of the physical subroutines 
which checks to make sure that the output of the subroutine 
is physically reasonable. See section 2.2.2.3. 


Used to show the name of the variable to which an integer 
pointer into COMMON /VAR/ corresponds. Used in conjunction 
with VARNAM. 

SUBROUTINE. 


Main controlling program of CFC VI. Calls the other sections 
of the program to control the integration and to print out the 
results. See section 4.1.1. 

PROGRAM in the control section of the interface. 


Maximum number of objects allowed in each room. 
INTEGER PARAMETER. 


The maximum integration order which the numerics is allowed 
to use. . 
INTEGER PARAMETER in MARK6. 


Maximum number of rooms allowed. 
INTEGER PARAMETER. 


The maximum number of variables which are allowed in the 
system. 
INTEGER PARAMETER 
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MAXVTS 


MAXWLS 


MNXCNC 


MNXPRD 


MO 


NETMAS 


NEWSTP 


non-eliminated 
variable 
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Maximum number of vents allowed. 
INTEGER PARAMETER. 


Maximum number of walls allowed. 
INTEGER PARAMETER. 


Calculates the mass concentration of carbon monoxide in a 
room. 
SUBROUTINE in the physics module. 


Calculates the rate at which the fire over an object produces 
carbon monoxide. : 
SUBROUTINE in the physics module. 


A synonymn for MAXOBS. 
INTEGER PARAMETER. 


A synonymn for MAXRMS. 
INTEGER PARAMETER. 


A synonym for MAXVTS. 
INTEGER PARAMETER. 


A synonymn for MAXWLS. 
INTEGER PARAMETER. 


Checks for the conservation of mass in a room. 
SUBROUTINE in the physics module. 


(NEW STeP.) Performs initializations in the interface at the 
beginning of each time step. See section 4.1.5. 
SUBROUTINE in control section of interface. 


Any variable which is not an eliminated variable. See section 
4.2.3.1. 
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non-functional sub- 


routine 


NROOMS 


NT 


numerics mod- 
ule 


NVARS 


NVENTS 


NWALLS 


NWSTAT 


NXTDSK 


NXTTY 


A physical subroutine whose output depends on its previous 
inputs as well as on its immediate inputs. See section 2.2.2.4. 


Number of rooms in the building. 
INTEGER in COMMON /BLDNG/. 


The number of the time steps since the start of the run. 
Local INTEGER in MARK6. 


The section of the program which is responsible for solving the 
system of algebraic and differential equations which describe 
the progress of a fire. 


Size of section 1 of COMMON /VAR/. 
INTEGER PARAMETER. 


Number of vents in the building. 
INTEGER in COMMON /BLDNG/. 


Number of walls in the building. 
INTEGER in COMMON /BLDNG/. 


(NeW STATe.) Reinitializes the system of equations whenever 
there is a change in the system of equations to be solved. See 
section 4.1.3. 

SUBROUTINE in control section of interface. 


The time at which the next output to the disk file FIRE.DAT 
is scheduled to occur. 
Local REAL¥8 in MARK6. 


The time at which the next output to the terminal screen is 
scheduled to occur. 
Local REAL*8 in MARK6. 
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OBJS (R) 


OBVARS 


OLDORD 


OLDTIM 


OLDX 


ORDER 


output variable 


OXYBRN 


OXYCNC 


PACKX 
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Number of objects in the room. 
INTEGER in COMMON /BLDNG/. 


Number of arrays in section 1 of /VAR/ indexed XXXXX (MAXOBS, 
MAXRMS). 
INTEGER PARAMETER. 


Temporarily stores the order of the integration method so that 
the program can restart with that method if the integrator 
should step beyond an unplanned discontinuity. 

Local INTEGER in MARK6. 


The time at which the previous time step values in common 
/BLDNG/ were taken. 
REAL*8 in COMMON /BLDNG/. 


Used to save the values in X at the beginning of each time step 
in case they need to be recovered for any reason before the 
completion of the time step. 

Local REAL*8 array in MARK6 and ALDIFF. 


The order of the integration method. 
Local INTEGER in ALDIFF and MARK6. 


In the context of a particular physical subroutine, a function 
value which is returned by that subroutine. See section 2.2.2.1. 


Calculates the mass rate of oxygen burning in the plume of an 
object. 
SUBROUTINE in the physics module. 


Calculates the mass concentration of oxygen in the hot layer 
of a room. 
SUBROUTINE in the physics module. 


Packs variables from section 1 of /VAR/ into the vector of 
variables which is sent to the numerics. 
SUBROUTINE in communications section of interface. 
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parameter Any input to a physical subroutine whose value cannot change 
in the middle of a time step. (Not to be confused with the 
PARAMETERs defined by the FORTRAN statement with that 
name.) Parameters are to be contrasted with variables whose 
value does change in the middle of a time step. See section 
orita: 


PCKPTR Tells subroutine FUNCT the order in which the results of a 
function call are to be packed from COMMON /VAR/ into the 
array Z before returning to the numerics. The integer values in 
this common block are set up by SETSYS each time the system 
of equations to be solved changes. 

COMMON block. 


PENTA Used by ABSRB3. 
SUBROUTINE in the physics module. 


PHI(0,R) Net flux incident on the surface of object 0 in room R. 
REAL*8 in section 1 of COMMON /VAR/. 


physical equation An equation which describes some physical process occurring 
in a fire. 


physical subrou- 
tine A FORTRAN SUBROUTINE which implements one or more 
physical equations. 


physical variable A variable for which the numerics must solve. 


physics module The collection of all the physical subroutines. See chapter 2. 


86 Appendix C: Dictionary 


planned discon- halk 
tinuity A discontinuity whose time of occurrence is known before the 
integration reaches that time. In other words, a discontinuity 
which the program can plan for. See section 4.1.1. 
planned interrup- 
tion An interruption whose time of occurrence is known before the 
integration reaches that time. In other words, an interruption 
which the program can plan for. See section 4.1.1. 
PLMHT Calculates the height of the plume over an object. 
SUBROUTINE in the physics module. 
POOL Used as a value for FIRTYP(0,R) to indicate that the object is 
a pool fire. 
INTEGER PARAMETER. 
POWL Calculates the rate of change of total energy of the hot layer. 
SUBROUTINE in the physics module. 
PPLR (R) Power from all plumes in room R to the hot layer in that room. 


REAL*8 in section 1 of COMMON /VAR/. 


previous time- 

step value The value of a physical variable at the end of the previous 
time step. Currently, a few of the physical subroutines take 
these values as input parameters. Previous time step values 
are stored in COMMON /BLDNG/. 


PSI0(0,R) Initial value for the semi-apex angle of the flame over the 
object. this value is used whenever the fire of the object is 
not oxygen starved. 

REAL*8 in COMMON /BLDNG/. 


PULRAD Calculates the radius of a pool fire. 
SUBROUTINE in the physics module. 


Appendix C: Dictionary 87 


PUTPTR 


PWRPLR 


PYRO 


PYROL 


QF (0,R) 


QVAP(0,R) 


RADFW 


RADLW 


RDNL 


RDNP 


RECAP 


Used repeatedly by SETSYS to perform a recurring pattern of 
assignments. 
SUBROUTINE. 


Calculates the power absorbed by the layer due to radiation 
from the plume of an object. 
SUBROUTINE in the physics module. 


Calculates the mass pyrolysis rate of a burning object. 
SUBROUTINE in the physics module. 


Used as a value for STATE(O,R) to indicate that the object is 


pyrolyzing. 
INTEGER PARAMETER. 


Heat of combustion as fuel of object 0 in room R. 
REAL*8 in COMMON /BLDNG/. 


Heat of vaporization of the object. 
REAL*8 in COMMON /BLDNG/. 


Calculates the radiative flux from the plume of an object to 
the walls in the same room. 
SUBROUTINE in the physics module. 


Calculates the radiative flux from the hot layer to the walls. 
SUBROUTINE in the physics module. 


Calculates the energy losses and gains of the hot layer due to 
radiation. 
SUBROUTINE in the physics module. 


Calculates the total radiative heat loss from the flame of an 
object. 
SUBROUTINE in the physics module. 


Writes a recapitulation of the input data to an output device. 
SUBROUTINE. 
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relevant variable A variable which is relevant to the fire currently being modelled. 
See section 4.1.4. 


RESULT Integer which gives information on the results of the just com- 
pleted call to the integrator. See section 3.1.3. 
Local INTEGER in ALDIFF and MARK6. 


RMVARS Number of arrays in section 1 of /VAR/ indexed XXXXX (MAXRMS). 
INTEGER PARAMETER. 


RNFFO2 Calculates the radiative flux absorbed by the flame of one 
object from the flame of another. 
SUBROUTINE in the physics module. 


RNLH Calculates the flux from the hot layer to the horizontal surface 
of an object. uy 
SUBROUTINE in the physics module. 


RNLV Calculates the flux from the hot layer to the vertical surface 
of an object. 
SUBROUTINE in the physics module. 


RNPOO1 Calculates the flux from the plume of one object to the surface 
of an object. 
SUBROUTINE in the physics module. 


RNWOO2 Calculates the flux from the walls in a room to an object in 
that room. 
SUBROUTINE in the physics module. 


root- finding equa- 
tion An equation in the system of the form given in (case a) of 
equation (2.5). 


root-finding vari- 
able A variable whose progress is described by a root-finding equa- 
tion. A component of the vector U in equation (1.1). 
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ROUND 


RTIME 


SCALE 


SCALER 


SCALER 


SCRICH 


SDVARS 


SETSYS 


- SIZE 


SMOLD 


SOLVE 


Smooths the singularity at 0 in the vent flow functions. 
REAL*8 FUNCTION. 


Returns the CPU time elapsed since the beginning of the run. 
(For use on Harvard VAX 11/780 only.) 
SUBROUTINE. 


Provides a magnitude of a physically significant value for each 
of the variables. 
SUBROUTINE in control section of interface. 


Common block corresponding to COMMON /VAR/ which contains 

one number for each variable representing the magnitude of 

the smallest physically significant value of that variable. 
COMMON block. 


Block data program section which initializes the values in COMMON 
/SCALER/. 
BLOCK DATA. 


Used by ABSRB3. 
SUBROUTINE in the physics module. 


Number of arrays in section 1 of /VAR/ indexed XXXxXx(2, 
MAXWLS). 
INTEGER PARAMETER. 


Selects the equations which are in the system during the various 
phases of the run. See section 4.1.4. 
SUBROUTINE in control section of interface. 


The number of equations in the system. 
Local INTEGER in interface and numerics. 


Integer parameter used as a value for STATE(O,R) to indicate 
that the object is smoldering. 
INTEGER PARAMETER. 


Solves a system of linear equations. 
SUBROUTINE in numerics module. 
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SOOT 


SPECFL 


SPREAD 


STATE(O,R) 


STRTIUP 


SUTCNC 


SUTPRD 


SYSPTR 


system 


system of equa- 
tions 
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Used by ABSRB3. 
SUBROUTINE in the physics module. 


Calculates the mass flow of each of the gaseous species through 
@ vent. 
SUBROUTINE in the physics module. 


Calculates the rate of spread of a growing fire. 
SUBROUTINE in the physics module. 


The burning state of object 0 in room R (e.g. flaming, heating, 
etc.) 
INTEGER in COMMON /BLDNG/. 


Value used by NWSTAT to indicate the start of a run. 
INTEGER PARAMETER. 


Calculates the mass concentration of soot in the hot layer of 
a room. 
SUBROUTINE in the physics module. 


Calculates the rate at which the fire on an object produces 
soot. 
SUBROUTINE in the physics module. 


Common block which tells subroutine FUNCT the order in which 
variables are stored in the linear array X which is sent from 
the numerics. See sections 4.1.4, 4.2.2, and the introduction 
to chapter 4. 

COMMON block. 


See system of equations. 


The set of n simultaneous equations in n variables which the 
numerics must solve in order to model a particular fire. 
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TCO2M 
TDISC 
TELZD (R) 


TELZR(R) 
TELZZ(R) 


TENPLM 
TEOSM (2) 
TEOZZ(0,R) 
ae 


TEPZR (0,R) 


Calculates the time derivative of the mass of carbon dioxide 
in the hot layer. 
SUBROUTINE in the physics module. 


When UNPINT detects an unplanned discontinuity, it sets this 
variable to the time of its occurrence. 
Local REAL*8 in UNPINT and MARK6. 


Rate of change of layer energy due to the convective heating 
of the walls. 
REAL*8 in section 1 of COMMON /VAR/. 


Total radiative power to the hot layer in room R. 
REAL*8 in section 1 of COMMON /VAR/. 


Time derivative of the energy of the hot layer in room R. 
REAL*8 in section 1 of COMMON /VAR/. 


Calculates the rate at which the plume over an object adds 
energy to the hot layer. 
SUBROUTINE in the physics module. 


Total rate at which all objects in room R release chemical 
energy. 
REAL*8 in section 1 of COMMON /VAR/. 


Negative of chemical energy release rate of object 0 in room 
R. 
REAL*8 in section 1 of COMMON /VAR/. 


Rate at which all the plumes in room R add energy to the hot 
layer in that room. 
REAL*6 in section 1 of COMMON /VAR/. 


Radiative power output from the plume of object 0 in room R. 
REAL*8 in section 1 of COMMON /VAR/. 
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TEUSM (R) 


TH20M 


TIME 


TINT 


TMDSM (R) 


TMDZF (R) 


TMFGZ(0,R) 


TMLZZ(R) 


TMMZF (R) 


TMNET (R) 
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Rate at which energy flows out of the hot layer of room R 
through all the vents in that room. 
REAL*8 in section 1 of COMMON /VAR/. 


Calculates the time derivative of the mass of water in the hot 
layer. 
SUBROUTINE in the physics module. 


The fire time which has elapsed since the start of the run. 
Local REAL*8 in ALDIFF and CORECT. 


The next time at which the integration must be interrupted to 

facilitate some action which requires converged values for the 

variables (e.g. output or a change in the system of equations.) 
Local REAL*8 in MARK6. 


Rate at which mass flows out of the cold layer of room R 
through all the vents in that room. 
REAL*8 in section 1 of COMMON /VAR/. 


Rate at which carbon dioxide flows out of room R through 
vents. 
REAL*8 in section 1 of COMMON /VAR/. 


Mass flow rate of the fuel for the object (if the object is a gas 
burner. 
REAL*8 in COMMON /BLDNG/. 


Time derivative of the mass of the hot layer in room R. 
REAL*8 in section 1 of COMMON /VAR/. 


Rate at which carbon monoxide flows out of room R through 
vents. 
REAL*8 in section 1 of COMMON /VAR/. 


Net mass balance in room R. TMNET is the output of the mass 
balance equation coded in subroutine NETMAS. To ensure the 
conservation of mass, it should be 0 for any solution of the 
system of equations. 

REAL*8 in section 2 of COMMON /VAR/. 
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TMNXM 


TMODS (R) 


TMOMS (R) 


TMOSM (R) 


TMOSS (R) 


TMOWS (R) 


TMOZB(0,R) 


TMOZF (R) 


TMOZZ (0 ,R) 


TMOZZP (0, R) 


TMPL 


| Calculates the time derivative of the mass of carbon monox- 


ide in the hot layer. 
SUBROUTINE in the physics module. 


Total mass rate of evolution of carbon dioxide from all objects 
in room R. ‘3 
REAL*8 in section 1 of COMMON /VAR/. 


Total mass rate of evolution of carbon monoxide from all 
objects in room R. 
REAL*8 in section 1 of COMMON /VAR/. 


Negative of sum of mass pyrolysis rates of all objects in room 
R. Shy 
REAL*8 in section 1 of COMMON /VAR/. 


Total mass rate of evolution of soot from all objects in room 
R. 
REAL*8 in section 1 of COMMON /VAR/. 


Total mass rate of evolution of water vapor from all objects in 
room R. ee 
REAL*8 in section 1 of COMMON /VAR/. 


Mass rate of burning of the fuel from object 0 in room R. 
REAL*8 in section 1 of COMMON /VAR/. 


Rate at which oxygen flows out of room R through vents. 
REAL*#8 in section 1 of COMMON /VAR/. 


Negative of mass pyrolysis rate of object 0 in room R. 
REAL*8 in section 1 of COMMON /VAR/. 


Negative of mass pyrolysis rate at the previous time step. 
REAL*8 in COMMON /BLDNG/. 


Calculates the temperature of the hot layer in a room. 
SUBROUTINE in the physics module. 
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TMPO2 Calculates the surface temperature of an object. 
SUBROUTINE in the physics module. 
TMPSM(R) Rate at which the plumes in room R add mass to the hot layer, 
REAL*8 in section 1 of COMMON /VAR/. 
TMPW1 Calculates the temperature of the wall surface in a room. 
SUBROUTINE in the physics module. 
TMPZZ(0,R) Rate at which the plume over object 0 in room R adds mass 
to the hot layer. 
REAL*8 in section 1 of COMMON /VAR/. 
TMSL Calculates the time derivative of the mass of the hot layer. 
SUBROUTINE in the physics module. 
TMSPLM Calculates the rate at which the plume over an object adds 
mass to the hot layer. | 
SUBROUTINE in the physics module. 
TMSZF (R) Rate at which soot flows out of room R through vents. 
REAL*8 in section 1 of COMMON /VAR/. 
TMUSM (R) Rate at which mass flows out of the hot layer of room R 
through all the vents in that room. 
REAL*8 in section 1 of COMMON /VAR/. 
TMWZF (R) Rate at which water vapor flows out of room R through vents. 
REAL*8 in section 1 of COMMON /VAR/. 
TOXYM Calculates the rate of change of mass of oxygen in the hot 
layer. 
SUBROUTINE in the physics module. 
TPSIO(O,R) Tangent of the semi-apex angle of the flames over the object. 


This value is used whenever the fire of the object is not oxygen 
starved. 
REAL*8 in COMMON /BLDNG/. 
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transition func- 
tion 


TSUTM 


TTYOUT 
UNPACK 


UNPDSC 


UNPINT 


unplanned dis- 
continuity 


unplanned inter- 
ruption 


A function which provides the initial values when a new system 
of equations becomes applicable. See section 4.1.3. 


Calculates the rate of change of mass of soot in the hot layer. 
SUBROUTINE in the physics module. 


The time between outputs to the terminal screen during the 
Tun. 
Local REAL#*8 in MARK6. 


Unpacks the vector of variables sent from the numerics into 
section 1 of /VAR/. See section 4.2.1. 
SUBROUTINE in communications section of interface. 


The former name of UNPINT. May still appear in a version of — 
the code which differs in minor details from the one described 
in this document. | 


Determines if any unplanned interruptions have occurred upon 
the completion of a time step. See section 4.1.2. 
. SUBROUTINE in control section of interface. 


A discontinuity whose time of occurrence is not known before 
the integration has reached that time. See section 4.1.1. 


An interruption whose time of occurrence is not known before 
the integration has reached that time. See section 4.1.1. 


96 


VAR 


variable 


variable elimina- 
tion 


VARNAM 


VARSIZ 


varying 
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Contains one position corresponding to each possible variable 
of the model as well as a position designated to receive the out- 
put of each of the functions corresponding to the differential 
and root-finding equations. Used whenever it is necessary to 
refer to the variables by their semantically meaningful sym- 
bolic names rather than simply as components of a vector. See 
figure 10 and section 4.2.3.2. 
COMMON block. 


Any input to a physical subroutine whose value can change in 
the middle of a time step. The numerics tries to find the values 
of the variables as a function of time. (Not to be confused 
with FORTRAN variables.) Variables are to be contrasted 
with parameters which do not change in the middle of a time 
step. 


A procedure for finding a calling order for the physical sub- 
routines which yields a relatively small number of independent 
variables in the system. See section 4.2.3.1. 


Returns the FORTRAN name of the variable which corres- 
ponds to a position in COMMON /VAR/. (Used in conjunction 
with LOOKUP.) 

CHARACTER*10 FUNCTION. 


The number of variables in all of COMMON /VAR/ (including 
sections 2 and 3). 
INTEGER PARAMETER. 


A physical variable is varying if its value is changing during a 
particular portion of the progress of a given fire. 
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VARYNG 


VERSUN (8) 


VLINK(S,V) 
VMAOUT 
VMAZZ (R) 
VMOZZ(0,R) 
VMWZZ CW) 


VNICVR 


VIVARS 


Contains one seanyey location for every variable in section 1 of 
/VAR/. The g‘*integer in this common block indicates whether 
the value of the g** variable in common /VAR/ is undefined, 
fixed at a constant but well defined value, or varying according 
to one of several possible equations. The values in this common © 
block are set up by SETSYS each time the system of equations 
changes. See sections 4.1.4, 4.2.3, and the introduction to 
chapter 4. : 
COMMON block. 


Integer array indicating which version of subroutine § is in 
use (only for those physical subroutines with more than one 
version.) Currently, VERSUN is large enough to handle five 
subroutines with multiple versions. 

INTEGER in COMMON /BLDNG/. 


The room on the side 8 of the vent V. Note that if VLINK(S,V) = 
0 then the side § of the vent V is on the outside of the apse 
INTEGER in COMMON /BLDNG/. 


Density of the air outside the building. 
REAL¥8 in COMMON /BLDNG/. 


Density of the air in the cold layer of the room. 
REAL*8 in COMMON /BLDNG/. 


Density of the object. 
REAL*8 in COMMON /BLDNG/. 


Density of the wall material of wall W. 
REAL*8 in COMMON /BLDNG/. 


Calculates the area of the vent covered by the hot layer. 
SUBROUTINE in the physics module. 


Number of arrays in section 1 of /VAR/ indexed XXXXX (MAXVTS). 
INTEGER PARAMETER. 
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WLINK(S , W) 


WLVARS 


WRIT 


WTRCNC 


WIRPRD 


WVLS (R) 


XGAMAS (0 ,R) 


XGAMMA (0 ,R) 


The room on the side § of the wall W. Note that if WLINK(8,W) = 
O then the side & of the wall W is on the outside of the building. 
INTEGER in COMMON /BLDNG/. 


Number of arrays in section 1 of /VAR/ indexed XXXXX (MAXWLS). 
INTEGER PARAMETER. 


Writes out values for the variables as the run progresses. 
SUBROUTINE. 


Calculates the mass concentration of water vapor in the hot 
layer. 
SUBROUTINE in the physics module. 


Calculates the rate at which the fire of an object produces 
water vapor. 
SUBROUTINE in the physics module. 


(Width of the Vents in the hot Layer, Summed. Sorry about 
that.) Total length of the segments at which the interface 
between the hot and cold layers in room R intersects the vents 
in the walls of that room. 

REAL*8 in section 1 of COMMON /VAR/. 


Array used to store values for the variables and their time 
derivatives. 
Local REAL*8 array in numerics and interface. 


Stoichiometric air/fuel mass ratio for the burning of the object. 
REAL*8 in COMMON /BLDNG/. 


Air/fuel mass ratio in a free burn of the object. 
REAL*8 in COMMON /BLDNG/. 


Array used to store a significant value for each variable in the 
system (against which relative error checking is performed). 
Local REAL*8 array in ALDIFF and MARK6. 
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YCcO2 


YO2 


ZBVZZ(V) 


ZC0ZZ(0,R) 


ZCWZZ (W) 


ZELZZ(R) 


ZERO 


ZG0ZZ (0 ,R) 


ZGWZZ (W) 


ZHFZZ(0,R) 


ZHFZZP (0,R) 


ZHLZZ (R) 


Mass concentration of carbon dioxide in the outside air. 
REAL*8 in COMMON /BLDNG/. 


Mass concentration of oxygen in the outside air. 
REAL*8 in COMMON /BLDNG/. 


Width of vent V. 
REAL*8 in COMMON /BLDNG/. 


Specific heat of the object. 
REAL*8 in COMMON /BLDNG/. 


Specific heat of the wall W. 
REAL*8 in COMMON /BLDNG/. 


Total energy of the hot layer in room R. 
REAL*8 in section 1 of COMMON /VAR/. 


One dimensional array used to hold the results of a function 
call which the numerics places through the communications 
section of the interface. 

Local REAL*8 array in ALDIFF. 


Thermal diffusivity of the object. 
REAL*8 in COMMON /BLDNG/. 


Thermal diffusivity of the wall W. 
REAL*8 in COMMON /BLDNG/. 


Height of the flames on object 0 in room R. 
REAL*8 in section 1 of COMMON /VAR/. 


Flame height of object at the previous time step. 
REAL*8 in COMMON /BLDNG/. 


Depth of the hot layer in room R, measured from the ceiling. 
REAL*8 in section 1 of COMMON /VAR/. 
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ZHOZZ(0,R) 
ZHPZZ(0,R) 
ZHRZZ(R) 


ZHTZZ(V) 


ZHVZZ.(V) 
2.5022(0,R) 
ZIWZZ.(W) 
ZKAOQUT 
ZKAZZ(R) 
ZKFZZ 
2KLZZ(R) 
2KOIG(0,R) 


ZKOPY (0, R) 
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The height of the object from the floor. 
REAL*¥8 in COMMON /BLDNG/. 


Height of the plume over object 0 in room R. 
REAL*8 in section 1 of COMMON /VAR/. 


Height of the room. 
REAL*8 in COMMON /BLDNG/. 


Distance from the transom of the vent to the ceiling in the 
same room. 
REAL*8 in COMMON /BLDNG/. 


Height of the vent V (from base to transom). 
REAL*8 in COMMON /BLDNG/. 


Thermal conductivity of the object. 
REAL*8 in COMMON /BLDNG/. 


Thermal conductivity of the wall W. 
REAL*8 in COMMON /BLDNG/. 


Temperature of the ambient air outside the building. 
REAL*8 in COMMON /BLDNG/. 


Temperature of the cold layer in the room. 
REAL*8 in COMMON /BLDNG/. 


Temperature of the flames of the objects. 
REAL*8 in COMMON /BLDNG/. 


Temperature of the hot layer in room R. 
REAL*8 in section 1 of COMMON /VAR/. 


Ignition temperature of the object. 
REAL*8 in COMMON /BLDNG/. 


Temperature of the onset of pyrolysis for the object. 
REAL*8 in COMMON /BLDNG/. 
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ZK0ZZ(0,R) 


ZKOZZP (0, R) 


ZKWZZ (S , W) 


ZLAMDA 


ZLRZX (R) 


ZLRZY (R) 


ZMDL (R) 


ZMLZZ (R) 


ZMML (R) 


ZMOL (R) 


ZMOZO (0,R) 


ZM0ZZ(0,R) 


ZMOZZP (0, R) 


Temperature of the surface of object 0 in room R. 
REAL*8 in section 1 of COMMON /VAR/. 


Temperature of the upper surface of the object at the previous 
time step. 
REAL*8 in COMMON /BLDNG/. 


Temperature of side § of wall W. 
REAL*8 in section 1 of COMMON /VAR/. 


1./absorption coefficient of the flames of the objects. 
REAL*8 in COMMON /BLDNG/. 


Length of the room along the x-dimension. 
REAL*8 in COMMON /BLDNG/. 


Length of the room along the y-dimension. 
REAL*8 in COMMON /BLDNG/. 


Total mass of carbon dioxide in the hot layer in room R. 
REAL*8 in section 1 of COMMON /VAR/. 


Total mass of the hot layer in room R. 
REAL*8 in section 1 of COMMON /VAR/. 


Total mass of carbon monoxide in the hot layer in room R. 
REAL*8 in section 1 of COMMON /VAR/. 


Total mass of oxygen in the hot layer in room R. 
REAL*8 in section 1 of COMMON /VAR/. 


Initial fuel mass of the object. 
REAL*8 in COMMON /BLDNG/. 


Remaining fuel mass of object 0 in room R. 
REAL*8 in section 1 of COMMON /VAR/. 


Mass of the object at the previous time step. 
REAL*8 in COMMON /BLDNG/. 
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ZMSL (R) 


ZMWL (R) 


2NOZZ(0,R) 


ZNWZZ (W) 


ZOAZZ 


ZOWZM 


ZOWZN 


ZPRZZ(R) 


ZRFZI(0,R) 


ZRFZM (0, R) 


ZRFZO (0, R) 


ZRFZZ(0,R) 


ZRFZZP (0 ,R) 
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Total mass of soot in the hot layer in room R. 
REAL#*8 in section 1 of COMMON /VAR/. 


Total mass of water vapor in the hot layer in room R. 
REAL*8 in section 1 of COMMON /VAR/. 


The thickness of the object. 
REAL*8 in COMMON /BLDNG/. 


Thickness of the wall W. 
REAL*8 in COMMON /BLDNG/. 


Heat transfer coefficient for all objects. 
REAL*8 in COMMON /BLDNG/. 


Maximum heat transfer coefficient for walls. 
REAL*8 in COMMON /BLDNG/. 


Minimum heat transfer coefficient for walls. 
REAL#*8 in COMMON /BLDNG/. 


Atmospheric pressure in room R (above ambient). 
REAL*8 in section 1 of COMMON /VAR/. 


The radius of the object. 
REAL*8 in COMMON /BLDNG/. 


Maximum burning radius. 
REAL*8 in COMMON /BLDNG/. 


Initial burning radius of the object. 
REAL*8 in COMMON /BLDNG/. 


Radius of the fire on object 0 in room R. 
REAL*8 in section 1 of COMMON /VAR/. 


Radius of the fire of the object at the previous time step. 
REAL*8 in COMMON /BLDNG/. 
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ZTIG(O,R) 


ZUF ZZ 


ZULOUT 


ZULZZ (R) 


ZX0ZZ (0,R) 


ZYLDZ(R) 


ZYLMZ (R) 


ZYLOZ(R) 


ZYLSZ(R) 


ZYLWZ(R) 


ZY0ZZ(0,R) 


Time of ignition of the object. 
REAL*8 in COMMON /BLDNG/. 


Absorptivity of the flames of the objects. 
REAL*8 in COMMON /BLDNG/. 


Absorptivity of the air outside the building. 
REAL*8 in COMMON /BLDNG/. 


Infra-red absorptivity of the hot layer in room R. 
REAL#8 in section 1 of COMMON /VAR/. 


The x-coordinate of the object. 
REAL*8 in COMMON /BLDNG/. 


Mass concentration of carbon dioxide in the hot layer of room 
R. 
REAL*8 in section 1 of COMMON /VAR/. 


Mass concentration of carbon monoxide in the hot layer of 
room R. 
REAL*8 in section 1 of COMMON /VAR/. 


Mass concentration of oxygen in the hot layer of room R. 
REAL*8 in section 1 of COMMON /VAR/. 


Mass concentration of soot in the hot layer of room R. 
REAL*8 in section 1 of COMMON /VAR/. 


Mass concentration of water vapor in the hot layer of room R. 
REAL*8 in section 1 of COMMON /VAR/. : 


The y-coordinate of the object. 
REAL*8 in COMMON /BLDNG/. 
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